diff --git a/CHANGELOG.md b/CHANGELOG.md index f7eca04f..0c3f2e0c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **`ChaisemartinDHaultfoeuille.by_path` + non-binary integer treatment** — `by_path=k` now accepts integer-coded discrete treatment (D in Z, e.g. ordinal `{0, 1, 2}`); path tuples become integer-state tuples like `(0, 2, 2, 2)`. The previous `NotImplementedError` gate at `chaisemartin_dhaultfoeuille.py:1870` is replaced by a `ValueError` for continuous D (e.g. `D=1.5`) at fit-time per the no-silent-failures contract — the existing `int(round(float(v)))` cast in `_enumerate_treatment_paths` is now defensive (no-op for integer-coded D). Validated against R `did_multiplegt_dyn(..., by_path)` for D in `{0, 1, 2}` via the new `multi_path_reversible_by_path_non_binary` golden-value scenario (78 switchers, 3 paths, single-baseline custom DGP, F_g >= 4): per-path point estimates match R bit-exactly (rtol ~1e-9 on event horizons; rtol+atol envelope for placebo near-zero values), per-path SE inherits the documented cross-path cohort-sharing deviation (~5% rtol observed; SE_RTOL=0.15 envelope). **Deviation from R for D >= 10:** R's `did_multiplegt_by_path` derives the per-path baseline via `path_index$baseline_XX <- substr(path_index$path, 1, 1)`, which captures only the first character of the comma-separated path string (e.g. for `path = "12,12,..."` it captures `"1"` instead of `"12"`); this mis-allocates R's per-path control-pool subset for D >= 10. Python's tuple-key matching is correct in this regime — the per-path point estimates we compute are correct; R's per-path subset for the same path is buggy. The shipped parity scenario stays in `D in {0, 1, 2}` to avoid the R bug. R-parity test at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathNonBinary`; cross-surface invariants regression-tested at `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathNonBinary`. +- **New `paths_of_interest` kwarg on `ChaisemartinDHaultfoeuille`** for user-specified treatment-path subsets, alternative to `by_path=k`'s top-k automatic ranking. Mutually exclusive with `by_path`; setting both raises `ValueError` at `__init__` and `set_params` time. Each path tuple must be a list/tuple of `int` of length `L_max + 1` (uniformity validated at `__init__`; length match against `L_max + 1` validated at fit-time); `bool` and `np.bool_` are explicitly rejected, `np.integer` accepted and canonicalized to Python `int` for tuple-key consistency. Duplicates emit a `UserWarning` and are deduplicated; paths not observed in the panel emit a `UserWarning` and are omitted from `path_effects`. Paths appear in `results.path_effects` in the user-specified order, modulo deduplication and unobserved-path filtering. Composes with non-binary D and all downstream `by_path` surfaces (bootstrap, per-path placebos, per-path joint sup-t bands, `controls`, `trends_linear`, `trends_nonparam`) — mechanical filter on observed paths via the same `_enumerate_treatment_paths` call site, no methodology change. **Python-only API extension; no R equivalent** — R's `did_multiplegt_dyn(..., by_path=k)` only accepts a positive int (top-k) or `-1` (all paths). The `by_path` precondition gate at `chaisemartin_dhaultfoeuille.py:1118` (drop_larger_lower / L_max / `heterogeneity` / `design2` / `honest_did` / `survey_design` mutex) and the 11 `self.by_path is not None` activation branches in `fit()` were rerouted to fire under either selector. Validation + behavior + cross-feature regressions at `tests/test_chaisemartin_dhaultfoeuille.py::TestPathsOfInterest`. - **HAD `practitioner_next_steps()` handler + `llms-full.txt` reference section** (Phase 5). Adds `_handle_had` and `_handle_had_event_study` to `diff_diff/practitioner.py::_HANDLERS`, routing both `HeterogeneousAdoptionDiDResults` (single-period) and `HeterogeneousAdoptionDiDEventStudyResults` (event-study) through HAD-specific Baker et al. (2025) step guidance: `did_had_pretest_workflow` (step 3 — paper Section 4.2 step-2 closure on the event-study path), an estimand-difference routing nudge to `ContinuousDiD` (step 4 — fires when the user wants per-dose ATT(d) / ACRT(d) curves rather than HAD's WAS estimand and has never-treated controls; framed around estimand difference, NOT around the existence of untreated units, since HAD remains valid with a small never-treated share per REGISTRY § HeterogeneousAdoptionDiD edge cases and explicitly retains never-treated units on the staggered event-study path per paper Appendix B.2 / `had.py:1325`), `results.bandwidth_diagnostics` inspection on continuous designs and simultaneous (sup-t) `cband_*` reading on weighted event-study fits (step 6), per-horizon WAS event-study disaggregation (step 7), and the explicit design-auto-detection / last-cohort-only-WAS framing (step 8). Symmetric pair: `_handle_continuous` gains a Step-4 nudge to `HeterogeneousAdoptionDiD` for ContinuousDiD users on no-untreated panels (this direction is correct because ContinuousDiD's identification requires never-treated controls). Extends `_check_nan_att` with an ndarray branch via lazy `numpy` import for HAD's per-horizon `att` array; uses `np.all(np.isnan(arr))` semantics so partial-NaN arrays (legitimate event-study output under degenerate horizon-specific designs) do not over-fire the warning. Scalar path is bit-exact preserved across all 12 untouched handlers. Adds full HAD section + `HeterogeneousAdoptionDiDResults` / `HeterogeneousAdoptionDiDEventStudyResults` blocks + `## HAD Pretests` index covering all 7 pretest entry points + Choosing-an-Estimator row to `diff_diff/guides/llms-full.txt` (the bundled-in-wheel agent reference); the documented constructor + `fit()` signatures match the real `HeterogeneousAdoptionDiD.__init__` / `.fit` API exactly (verified by `inspect.signature`-based regression tests). Tightens the existing `Continuous treatment intensity` Choosing row to surface ATT(d) vs WAS as the estimand differentiator. `docs/doc-deps.yaml` updated to remove the `llms-full.txt` deferral note on `had.py` and add `llms-full.txt` entries to `had.py`, `had_pretests.py`, and `practitioner.py` blocks. Patch-level (additive on stable surfaces). 26 new tests (16 in `tests/test_practitioner.py::TestHADDispatch` + 9 in `tests/test_guides.py::TestLLMsFullHADCoverage` + 1 fixture-minimality regression locking the "handlers are STRING-ONLY at runtime" stability invariant). Closes the Phase 5 "agent surfaces" gap; T21 pretest tutorial and T22 weighted/survey tutorial remain queued as separate notebook PRs. ## [3.3.2] - 2026-04-26 diff --git a/benchmarks/R/generate_dcdh_dynr_test_values.R b/benchmarks/R/generate_dcdh_dynr_test_values.R index 1c237cd1..e6119720 100644 --- a/benchmarks/R/generate_dcdh_dynr_test_values.R +++ b/benchmarks/R/generate_dcdh_dynr_test_values.R @@ -927,6 +927,97 @@ scenarios$multi_path_reversible_by_path_trends_nonparam <- list( results = extract_dcdh_by_path(res18, n_effects = 3, n_placebos = 1) ) +# Scenario 19: by_path + non-binary integer treatment (D in {0, 1, 2}). +# Phase 3 Wave 3 #8 lift. Custom inline DGP (mirror Scenario 17 structure) +# with 3 single-baseline non-binary paths: low-dose sustained +# (0, 1, 1, 1), high-dose sustained (0, 2, 2, 2), and ramp-up +# (0, 1, 2, 2). All F_g >= 4 (defensive: avoids any pre-window boundary +# edge cases under future trends_lin combinations and matches Scenario 17). +# 78 switchers + 20 never-treated (D=0) + 20 always-treated (D=2) controls. +# n_periods=13, L_max=3. +# +# R's substr(path, 1, 1) baseline-derivation in did_multiplegt_by_path +# is correct for D in {0..9} (single-digit decimal); we stay in {0, 1, 2} +# so no R bug interferes. Python's tuple-key matching is correct +# regardless of D range. +cat(" Scenario 19: multi_path_reversible_by_path_non_binary\n") +{ + set.seed(119) + n_periods19 <- 13 + L_max19 <- 3 + target_paths19 <- list( + c(0L, 1L, 1L, 1L), # path 1, low-dose sustained (rank 1) + c(0L, 2L, 2L, 2L), # path 2, high-dose sustained (rank 2) + c(0L, 1L, 2L, 2L) # path 3, ramp-up (rank 3) + ) + fg_path_counts19 <- list( + list(F_g = 4L, path_idx = 1L, count = 18L), + list(F_g = 5L, path_idx = 1L, count = 14L), + list(F_g = 6L, path_idx = 2L, count = 14L), + list(F_g = 7L, path_idx = 2L, count = 12L), + list(F_g = 8L, path_idx = 3L, count = 12L), + list(F_g = 9L, path_idx = 3L, count = 8L) + ) + n_switchers19 <- sum(sapply(fg_path_counts19, function(x) x$count)) + stopifnot(n_switchers19 == 78L) + D19 <- matrix(0L, nrow = n_switchers19, ncol = n_periods19) + g19 <- 1L + for (entry in fg_path_counts19) { + F_g <- entry$F_g + target <- target_paths19[[entry$path_idx]] + n_here <- entry$count + for (k in seq_len(n_here)) { + if (F_g >= 3L) D19[g19, 1:(F_g - 2L)] <- 0L + for (j in 0:L_max19) D19[g19, F_g - 1L + j] <- target[j + 1L] + if (F_g + L_max19 <= n_periods19) { + D19[g19, (F_g + L_max19):n_periods19] <- target[L_max19 + 1L] + } + g19 <- g19 + 1L + } + } + # Append 20 never-treated (D=0) and 20 always-treated (D=2) controls + D19 <- rbind( + D19, + matrix(0L, nrow = 20L, ncol = n_periods19), + matrix(2L, nrow = 20L, ncol = n_periods19) + ) + n_total19 <- nrow(D19) + set.seed(119L) + group_fe19 <- rnorm(n_total19, 0, 2.0) + noise19 <- matrix(rnorm(n_total19 * n_periods19, 0, 0.5), + nrow = n_total19, ncol = n_periods19) + period_arr19 <- 0:(n_periods19 - 1L) + Y19 <- 10.0 + + matrix(group_fe19, nrow = n_total19, ncol = n_periods19) + + matrix(0.1 * period_arr19, nrow = n_total19, ncol = n_periods19, byrow = TRUE) + + 1.5 * D19 + + noise19 + d19 <- data.frame( + group = rep(seq_len(n_total19) - 1L, each = n_periods19), + period = rep(period_arr19, n_total19), + treatment = as.vector(t(D19)), + outcome = as.vector(t(Y19)) + ) + res19 <- did_multiplegt_dyn( + df = d19, outcome = "outcome", group = "group", time = "period", + treatment = "treatment", effects = 3, placebo = 1, by_path = 3, + ci_level = 95 + ) + scenarios$multi_path_reversible_by_path_non_binary <- list( + data = list( + group = as.numeric(d19$group), + period = as.numeric(d19$period), + treatment = as.numeric(d19$treatment), + outcome = as.numeric(d19$outcome) + ), + params = list(pattern = "single_baseline_multi_path_non_binary", + n_switcher_groups = 78L, n_realized_groups = 118L, + n_periods = 13L, seed = 119L, effects = 3, placebo = 1, + by_path = 3, ci_level = 95), + results = extract_dcdh_by_path(res19, n_effects = 3, n_placebos = 1) + ) +} + # --------------------------------------------------------------------------- # Write output # --------------------------------------------------------------------------- diff --git a/benchmarks/data/dcdh_dynr_golden_values.json b/benchmarks/data/dcdh_dynr_golden_values.json index ccc0a3d5..38b6a847 100644 --- a/benchmarks/data/dcdh_dynr_golden_values.json +++ b/benchmarks/data/dcdh_dynr_golden_values.json @@ -1300,6 +1300,143 @@ } ] } + }, + "multi_path_reversible_by_path_non_binary": { + "data": { + "group": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 71, 71, 71, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 76, 76, 76, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 79, 79, 79, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 82, 82, 82, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 83, 83, 83, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 86, 86, 86, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 102, 102, 102, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 103, 103, 103, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 104, 104, 104, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 107, 107, 107, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 108, 108, 108, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 109, 109, 109, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 113, 113, 113, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 114, 114, 114, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 117, 117, 117, 117, 117, 117, 117, 117], + "period": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], + "treatment": [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2], + "outcome": [5.5529442619, 6.156927173, 6.6707706304, 8.1432916815, 7.4323516067, 7.7140956247, 7.7632946841, 8.6254763535, 8.3597776755, 8.223141545, 8.180318544, 8.0041503186, 8.0021883778, 10.4956752712, 10.7008327224, 10.6320129542, 12.6213246145, 12.9000685162, 12.3408443357, 13.0557366696, 12.7120697282, 13.9273768144, 13.4513509537, 12.6328678723, 13.4705249374, 12.9117044684, 8.2531199663, 8.8711169095, 8.8940337746, 10.1788167709, 11.2666200793, 10.3184657214, 11.5285926565, 11.0822608472, 10.5368842007, 11.8785487447, 10.8852012079, 11.2475207018, 11.0043238688, 11.2683601183, 11.0513090077, 11.4988901391, 11.9001336355, 11.7835830811, 12.7595040156, 13.1380168069, 14.2884086945, 13.0885670645, 13.1345511082, 13.6076545078, 14.2953201489, 13.3435573873, 10.6464660953, 9.604989008, 10.9421373091, 12.5898132605, 11.6948303507, 12.2474700574, 12.2260140744, 12.131071744, 12.3093952221, 13.3101751362, 12.9865433189, 12.9163073698, 12.9299944361, 8.9154119787, 8.6799465242, 9.6921382553, 10.9308164058, 10.7893232463, 10.9921183855, 10.5611404307, 11.1261451104, 10.8001842313, 11.8236535629, 11.3242426209, 11.5273715252, 11.3783509903, 9.9269414184, 11.4258848645, 10.9783846177, 13.0942526694, 11.9222620785, 12.0334155522, 11.3883259813, 12.6645720294, 12.7581285368, 13.7059231734, 12.831935872, 13.1738684534, 13.8963321411, 5.7852838816, 4.7877667851, 5.7090419435, 7.8435764938, 6.8940074487, 6.5028708489, 7.1520632823, 6.9262008808, 7.8970978052, 7.812700529, 7.0710356808, 7.609998536, 8.2759983388, 10.8148440406, 12.3200515305, 11.9161718993, 14.0868310329, 13.5873791949, 14.2396251107, 14.2103750535, 14.3950996926, 14.2926315191, 14.0548255048, 13.3869407738, 14.6169860447, 13.8398661868, 10.8511942927, 10.8749883675, 10.4819846011, 12.8290793912, 12.7711038098, 12.2490728205, 12.4701860332, 13.8012762153, 12.9463844432, 12.9557086977, 12.2271809889, 12.9038886296, 13.13726656, 13.0547014058, 13.8116092102, 14.3273968085, 15.4198806316, 15.3012314147, 15.3775524984, 16.2716139253, 16.0157075856, 16.4927917098, 15.5790437569, 15.8887043617, 16.4282332904, 16.0546394736, 7.6515463806, 8.441454703, 8.2208259541, 10.8151890957, 9.9193404416, 10.533136181, 10.2465382314, 9.4461499127, 10.3639871557, 10.5562203626, 10.9732500075, 11.512706594, 10.892784061, 10.8583323628, 11.7324354206, 11.9001594377, 12.7355328111, 12.6916528065, 13.0365591956, 13.0492452969, 13.8230033292, 13.934161237, 12.5015998987, 13.3451448856, 13.721214858, 14.4362284246, 7.9813377739, 9.5066399723, 8.5677302664, 10.4432456562, 10.4184189859, 10.1733642886, 9.8976084566, 10.6853056941, 11.5113721858, 11.0176445211, 11.3152098433, 12.1803455391, 12.3452567074, 7.6331450185, 8.2093247233, 8.4818386217, 10.6337713245, 9.8877167344, 9.493966436, 10.5436141766, 10.1060178419, 10.2386325993, 10.4555261704, 10.9559109084, 10.847196989, 10.513027499, 10.9567992108, 10.4638478834, 10.4450336796, 11.4216447264, 11.7393450072, 12.9024451291, 12.2621370958, 12.5791745655, 12.7945249308, 12.3240998479, 13.0779725599, 13.4044651956, 11.8447578016, 8.3228607619, 8.9038035375, 10.0466507673, 11.1881600237, 11.0811023012, 10.9087363528, 11.1908306303, 10.7125072026, 11.0424888987, 11.9509489851, 10.7202516556, 11.1707456744, 11.6863681776, 10.2900125356, 9.9942576394, 10.670886949, 12.5205400141, 12.2760892014, 11.6362925367, 12.3264641363, 13.1617046939, 13.0913822179, 12.3517721688, 13.0201750068, 13.2391551791, 13.2316777877, 4.1382995464, 3.8027044415, 5.0056889214, 4.2965142646, 5.4938366709, 4.9549193906, 5.1736561835, 6.9599742667, 5.5202824919, 5.9695804004, 6.5464030749, 5.4313861438, 5.6619650332, 12.0162745546, 12.9102093713, 13.9885272852, 12.6989023744, 15.8232350854, 13.7291125545, 14.6237425912, 15.0114832463, 14.6181743554, 14.6649179396, 15.0803252233, 15.3625329787, 15.6406160247, 11.8892720911, 12.0979372462, 11.4736429724, 12.014214575, 14.302178722, 13.9022894351, 13.5824717688, 13.8723011082, 14.2124929485, 14.3086070779, 15.2029465707, 13.725152849, 14.1304087158, 8.8634223075, 8.9909593051, 9.464655569, 9.4521589979, 10.2182058649, 11.0579496468, 11.3164298108, 11.1953772137, 11.1568627018, 11.8090537615, 11.6261925899, 11.129529477, 11.9836078661, 10.3526463825, 9.5692371366, 9.790125867, 10.2965334593, 11.3866037267, 11.2516967067, 11.4398346606, 11.5975899265, 11.6035440405, 12.0403055423, 12.1525280735, 11.7598334329, 12.7797906589, 6.1368993731, 6.3214414418, 5.8359049884, 6.2427029732, 7.982466959, 7.9666707187, 8.2640361619, 8.377783177, 7.5809864822, 7.7128677077, 7.8164677209, 8.5821178194, 7.027058568, 10.8781089188, 10.3006595359, 10.435311488, 11.3405020612, 11.0865446198, 11.9944511893, 11.5024171643, 12.5236053722, 12.285306064, 13.2055327423, 12.3647076707, 12.8980139269, 14.2250133172, 9.1345732637, 9.4063058338, 10.2402174305, 9.8418300361, 10.7088053172, 11.324818406, 11.8709625776, 11.3909176994, 12.1947605163, 12.4059040821, 11.6142703549, 12.9367852792, 12.9601708608, 7.2743263233, 7.6206127983, 8.3680408686, 7.8841166705, 9.7613910744, 10.4370019227, 10.5457648913, 10.3932252796, 10.4813509287, 9.71414913, 10.4023135358, 10.5753440572, 9.7811949964, 10.658840068, 11.4638413736, 11.4524777624, 12.3656308817, 13.1734565108, 12.9507259761, 13.1758442972, 14.1117107849, 13.7832373746, 13.417214655, 13.0631548284, 14.3453924348, 14.5592368427, 11.2454083509, 10.6944270147, 11.4046854784, 10.4196622695, 13.124578303, 13.2996293119, 12.8530943874, 12.8138176723, 12.9721338747, 13.1104511641, 14.9414891519, 13.4154745859, 13.9699866221, 10.4278647804, 10.7159720588, 10.2908785626, 10.3208730937, 12.7109846263, 12.1957500861, 12.037099632, 12.1035014287, 12.10919479, 12.4770405823, 12.9003538204, 13.2798146513, 13.0334480937, 10.6606790521, 9.4629759899, 9.9226420793, 10.6755396203, 12.0809255721, 13.0681109036, 11.9331392989, 13.2176031468, 13.0900218855, 12.6934866682, 12.375648565, 13.2993595746, 13.3478618121, 7.7187791924, 7.5361802936, 8.4475599477, 9.445681525, 9.7627734638, 10.5197373304, 9.6399867055, 10.3090292604, 10.0367071983, 9.888588554, 10.9738800856, 10.2607701898, 10.5130131136, 11.3904665217, 11.2217235951, 10.0849064751, 10.2867422805, 11.2325128054, 13.8824354497, 13.7255790101, 13.8133678629, 14.9120264004, 15.3454802504, 15.1145266802, 15.0891199285, 13.9729457173, 12.083720194, 12.2404851901, 12.3044859955, 12.3056885555, 13.1843699838, 15.3460757975, 15.5007552184, 15.4842931277, 16.1656120757, 15.4832255218, 16.4521636263, 15.777030394, 16.5466736799, 6.4129469636, 6.2958740534, 6.901104925, 6.924374479, 7.2407848213, 10.8987502644, 11.0584999789, 10.9310870349, 10.5645156131, 10.7228860109, 11.0107298792, 11.0735871204, 11.2610323385, 11.7524630742, 10.5389532507, 12.0789420768, 12.7921734529, 12.1825044948, 14.8669310011, 14.6252162924, 15.0064162597, 14.5581977991, 15.4328761498, 16.0226962325, 15.2538991841, 15.2156448088, 11.0011811813, 10.540955475, 10.7123535915, 11.0371519236, 11.2513900297, 14.1059653751, 13.3414886631, 14.369394526, 15.0928688465, 14.1964769505, 13.8924592398, 15.3615132616, 14.7812512329, 8.9442174355, 8.7353758991, 10.0215514873, 8.998843737, 9.6803540618, 12.4663026953, 12.1683827394, 11.7994070426, 12.2592415601, 13.3841802008, 13.7694226553, 12.2125494298, 12.7519493449, 11.6973944685, 11.612613791, 11.4236926165, 12.5816335195, 12.6126449958, 14.9798903746, 14.7686113463, 14.8735740328, 14.8206340166, 15.2069148481, 15.1416179314, 15.976877366, 16.1833039198, 10.205268761, 10.3813462545, 11.1135981147, 11.070799731, 11.5893568809, 14.13371698, 14.263237822, 14.7521100159, 15.3074442584, 15.1596967353, 15.4467091833, 15.9837015951, 15.0367564647, 11.7149029254, 10.8034907224, 10.6504609168, 11.2689403059, 11.812603087, 14.6416090075, 14.0672973637, 14.8974158752, 15.9285321398, 15.2610964335, 14.9736753301, 14.7753756457, 15.4623999009, 9.3628665578, 9.7234752706, 9.7828291242, 9.7379680732, 10.0889300468, 11.6144181289, 13.102643687, 12.9998425433, 13.2752516406, 13.8265724829, 12.5966105864, 14.0232641932, 12.9678349699, 12.9162510994, 12.6340894708, 14.1725265636, 13.700431448, 13.9020097389, 16.1811877383, 17.0532926566, 16.7460693518, 16.6233951838, 16.8979967385, 16.9244730814, 16.9445487702, 17.5390650612, 8.8946934836, 9.4498980886, 8.7082451488, 9.3695307253, 9.6352320328, 11.817349757, 11.8956000845, 12.4509657889, 12.5105938136, 12.137238643, 11.8100689742, 12.3819750286, 12.7467956817, 10.1441525602, 10.380254364, 11.0606131928, 10.4810263222, 11.620677904, 13.5926966194, 14.3528196382, 13.0364163884, 13.9981610938, 13.5776360772, 15.2632830916, 15.8792780565, 15.0485159799, 13.7957025431, 14.6655979042, 14.8558351974, 15.244837639, 14.1668638806, 17.9033209722, 18.1443726897, 17.7799956555, 17.074591048, 18.0476585757, 19.2587033057, 18.5005584428, 18.8879024564, 10.9650229779, 10.8314154581, 9.7214569016, 10.8358876567, 11.44338465, 10.9903036071, 13.9455194378, 14.3525999359, 14.2042495774, 14.7832489612, 14.3472736458, 14.7736315587, 15.0392783219, 7.7186419121, 9.2098118187, 7.9932691868, 8.6683945259, 8.8726750455, 8.9172281038, 11.4939011262, 12.1950320993, 11.3960972834, 13.0098145584, 12.9716986269, 12.4596330208, 11.9994017797, 10.6376739845, 11.0130987999, 11.312090153, 11.6615479403, 11.8428602685, 11.431940196, 13.6795925963, 14.3788361067, 15.2378433189, 13.6468806302, 14.7906525924, 15.3049737605, 15.476610981, 5.7035340214, 6.7004400299, 7.4194813663, 6.7982788605, 7.9811707338, 7.081903396, 10.3469837894, 11.320666984, 10.9731845463, 10.8897300461, 10.5447259203, 11.232975723, 11.2042337895, 6.4106618855, 6.5906703158, 7.4390826456, 7.1218774648, 8.1009647218, 6.84804322, 10.3321529534, 9.1748997961, 10.5053769811, 10.7415287436, 10.2366068725, 10.6646781245, 10.4999794534, 11.9971095563, 12.2177062556, 11.6623938658, 12.7575771901, 12.5014382965, 12.1760429564, 15.13713986, 15.0871198194, 15.0813294232, 15.4804938411, 15.9265722061, 16.2463245923, 16.2234777739, 8.3321532831, 9.2882089502, 9.1488746454, 10.0681870418, 8.5302443104, 9.1339048269, 12.4489300243, 12.4877757183, 12.7145511679, 13.1458299845, 12.6183015274, 12.8371178753, 13.3920850692, 10.528583888, 10.4414682044, 10.6359836552, 11.5485615586, 10.9823175773, 11.4032932215, 14.8176983764, 14.163005259, 13.5682191036, 14.9763651786, 14.3852993219, 14.3141363119, 14.2743613573, 9.9070551741, 9.8597446086, 9.0616193724, 9.9606537335, 9.3017847839, 9.5854568091, 12.6711661282, 13.3937357591, 12.1532485281, 12.9803160491, 12.7169751642, 14.5339176447, 13.0697492346, 7.0696013944, 8.5035252391, 8.2986451035, 8.0612102886, 7.8275416795, 8.0243967112, 11.4122171186, 11.1959650594, 11.8125743615, 11.9422735884, 12.0578014253, 13.1616746628, 11.8179924386, 9.6434574371, 10.4873090751, 9.6630404886, 9.85015153, 10.2290097308, 10.412859459, 13.1397968508, 13.5730270739, 13.58038109, 13.62984411, 14.0397101195, 13.0617692895, 13.5492809142, 10.3561316452, 10.1861110751, 10.1892449871, 10.2685083685, 10.3625791342, 10.8391026998, 14.3846657179, 14.7538372895, 14.7299097121, 14.2030557585, 14.43652827, 13.9691059086, 14.3631771625, 10.0509937821, 9.1976340328, 10.5064851881, 11.3560502582, 10.3958512206, 10.8242875264, 10.199508338, 11.0721783444, 14.2011837605, 14.9460836538, 13.5595864735, 13.7068133659, 13.7798031468, 6.4578381007, 6.2841933968, 6.4766451579, 7.126698767, 6.6194569904, 6.4329305959, 7.436586626, 8.2030055073, 9.5629734446, 9.8176945176, 10.2677309882, 9.6178316797, 10.5565141928, 11.2019550173, 11.1099262211, 10.086956578, 10.4509110072, 11.1528998687, 10.9881937556, 11.3835654073, 13.3752237068, 14.4289307836, 14.955103696, 14.6142469099, 14.4328087077, 15.0313201528, 10.9593762192, 11.7010676256, 11.3537173902, 11.8390482919, 11.3136354904, 12.3312774033, 11.6759968222, 13.5432182502, 15.6416776251, 15.7724851162, 15.7446223948, 16.2252630099, 15.1638476824, 12.6032190012, 12.7080336248, 12.6871928878, 12.9161923969, 12.6450011982, 12.5486742937, 12.047372472, 14.1085936652, 16.4673118597, 15.8797526939, 15.4441630533, 17.1175166166, 16.4078387472, 6.44433104, 6.6046795453, 7.6857008347, 7.0349174211, 8.2615279946, 8.079487139, 7.7719218661, 9.0015294317, 10.7898324695, 11.5633506728, 10.8129376913, 12.0848551408, 11.6976973571, 8.0136928692, 7.6808573108, 8.2972630961, 8.6565638305, 8.5443040376, 8.6411490615, 8.735142359, 10.4987545944, 11.3897814302, 12.5049704283, 11.7864061037, 12.5706186403, 12.1464766952, 10.9150668709, 11.1683255699, 10.5167607057, 11.0262485691, 11.2663379937, 10.584672605, 11.1078221707, 12.7230626412, 13.9870967477, 14.2086174289, 14.4287988528, 13.7996893426, 14.5217528749, 12.262573502, 11.925454021, 11.3352927929, 11.3470512808, 11.0823171854, 11.6372893285, 11.9854190602, 13.2331844763, 15.2481243272, 15.1582631184, 15.8431494021, 15.1239757151, 15.519406832, 12.112736165, 12.0932629664, 10.4995798152, 11.5312956384, 11.8951530053, 11.5181887313, 12.5862324524, 13.4809846125, 14.2906942806, 15.4845326193, 15.2775372503, 15.9906608747, 15.5123185611, 10.4324150208, 10.9669452615, 9.2004843606, 9.9196252728, 10.9244317478, 10.3134388781, 9.7721436209, 11.6255896018, 13.8890833975, 13.4255932247, 14.2702484378, 13.2954697747, 14.0698068422, 6.2549390684, 6.3285172336, 5.9094888466, 6.2497626796, 5.7273267029, 6.2668030449, 6.6726840996, 8.29021421, 10.8317229448, 9.4617771496, 10.2015834309, 10.3594292127, 9.9206246465, 3.8729059276, 5.1032063749, 4.7393154973, 4.965059281, 4.9112950564, 5.6736853961, 4.8851157199, 5.9085041257, 6.1033168011, 7.9320134511, 8.0710776681, 8.4733137135, 8.6634503445, 8.2333211946, 8.8849395662, 8.4158537717, 8.8195667305, 8.7141105792, 9.4841683481, 9.6937438523, 9.5768077083, 11.1956477494, 12.4597477442, 14.0604098141, 12.8363162584, 13.4096328605, 15.4510526192, 15.6338600508, 15.3439614, 15.545247528, 15.9060931871, 17.090082964, 15.9327188571, 16.6342892453, 18.4606312504, 19.8194333515, 19.148061831, 18.5238831851, 20.1252006052, 9.2597746261, 9.2825371557, 10.264921785, 9.9332966016, 9.2426929428, 10.7887917379, 10.2008482315, 11.7274562906, 12.202206483, 13.2799172503, 13.4640628398, 13.9828291886, 14.0777876383, 9.0885811883, 8.1935681046, 8.6741169445, 7.8861155803, 9.324512564, 7.9499583884, 8.6581382968, 8.843067005, 10.9177911561, 12.349667725, 11.6645330358, 12.0196318979, 12.7363160269, 10.2372852111, 11.2022097172, 10.6596720728, 10.3819594364, 10.676740047, 11.0009869083, 10.9404041971, 10.640478961, 12.9431393123, 13.9882584004, 13.9899278795, 14.8941962741, 14.575744495, 8.1494230797, 7.2252547835, 8.3176618619, 8.6156097984, 7.6799671312, 7.8269984016, 9.0341935269, 9.4877494253, 10.3573198745, 11.7538485622, 11.5329362939, 11.181252374, 12.200254274, 10.4414109857, 12.1114068989, 11.5445699364, 10.8289603465, 10.6329446802, 11.1526809258, 11.5085557871, 11.4307026286, 14.2718864359, 15.2505599401, 14.5888514164, 14.8052861201, 15.37933499, 10.6455407561, 10.7698971424, 11.3323585787, 11.6028342883, 10.8562874316, 11.2592151327, 10.9003328626, 10.2468885162, 10.8923342749, 11.493455065, 11.7201521496, 11.5530773296, 10.9532003832, 10.77932875, 11.8838693172, 11.357382095, 12.0396415044, 11.0880593493, 11.5105412552, 11.5559675904, 12.0178102281, 12.2101219169, 12.5560939297, 13.2448000009, 12.3509912808, 12.1460386388, 11.7885901305, 11.3908976578, 11.398424302, 11.5520333921, 11.0170241348, 12.2045348369, 11.159269782, 11.7019289777, 12.3664064572, 11.3194478071, 12.244961505, 11.8110465004, 12.5516960289, 12.6627092328, 12.9790314509, 12.9502488111, 12.5180228475, 12.6144397463, 12.7519397889, 12.3328231175, 12.8417614762, 13.3866685442, 12.8867128853, 13.3756456866, 14.1275523522, 12.4738075817, 11.6900601864, 11.9436697354, 12.4410597788, 10.3964085622, 12.1182113287, 11.9378623292, 12.1623002285, 12.2233516292, 12.372684249, 12.0500714383, 12.5088815428, 11.4543421647, 12.1440284933, 10.9515785066, 10.0662783572, 10.8951942727, 11.1313159736, 9.9812723136, 10.9248562731, 9.8533566907, 11.0909149394, 10.837619303, 10.5116157329, 11.1827362363, 10.9820671589, 11.3900195787, 7.79839664, 8.0252280564, 8.327815514, 7.1953311346, 7.3543738836, 8.046208592, 8.303950499, 8.9473088124, 8.1809799713, 8.5922869189, 8.6720426142, 8.7425623583, 8.246028714, 7.1755862544, 7.5913600858, 7.496272671, 8.3374518533, 7.8417669757, 8.0725415597, 7.8088316278, 8.5189557197, 6.7671073316, 8.0987086636, 8.1915897238, 8.018867779, 8.2485443274, 13.6818013238, 13.8007993357, 13.2692029041, 13.4551713121, 15.1237277246, 14.0647827788, 13.8687579892, 14.9313855656, 15.0510687997, 13.6533606668, 14.5958142369, 15.4177327021, 15.1932180404, 11.8287511378, 11.9568890933, 12.9939767778, 11.5751366337, 11.8750180382, 12.4649438414, 11.8816498014, 13.1214602404, 13.1176727944, 13.4800583631, 13.4913702387, 13.070061512, 12.9539471445, 10.3716136748, 10.2083263255, 9.7832463227, 10.3254005175, 10.8338157103, 11.2766684855, 10.7838076613, 10.5667018297, 11.3229045072, 9.6293960656, 11.5031219743, 10.4840972053, 11.3596752737, 7.5766526642, 7.9112211821, 7.2484664337, 7.0250735254, 7.5287636255, 8.6824218636, 8.8555344068, 7.9018599442, 8.192427109, 8.9909105614, 8.5883106146, 9.0585233917, 9.1951756463, 11.1539038322, 9.9020817372, 11.1304672348, 10.9078622834, 11.0022451771, 10.3264812362, 10.9096482242, 10.7138898281, 10.8341289912, 11.2796860954, 11.1310292388, 10.9041331068, 11.136487348, 10.1564461175, 9.3862887445, 10.1062405039, 9.8217075527, 10.9392843979, 10.6596731967, 10.939473531, 10.9922047009, 11.6028802155, 11.2918958164, 10.595375472, 11.6460114859, 10.6744884465, 10.0625294975, 10.0766393742, 9.7114369396, 10.6991415774, 10.4835980822, 10.7186547375, 10.3823253751, 10.5503857363, 10.5108304809, 11.5329683774, 10.8956816753, 11.343054347, 11.9023066354, 8.2966081169, 7.8767677741, 7.4955157472, 8.0916884965, 8.2025139581, 8.2380024126, 7.9519594603, 9.3252826258, 8.3330664871, 8.6208058872, 8.9071824273, 9.4799526775, 8.7808248403, 5.8185719811, 6.0098168518, 5.1598380542, 6.048613675, 6.2219264182, 5.6977113836, 6.0196039675, 6.7293741174, 5.8439839673, 6.727503313, 6.8456750489, 7.2418903609, 6.5159804324, 8.357812333, 9.2810336306, 8.6787696068, 8.1640309289, 8.0516837827, 8.384930169, 8.7787664101, 9.9823057828, 10.0816273435, 8.8587326188, 9.201553707, 9.9421723188, 9.8093913109, 6.9673205886, 7.1839694534, 7.3829792703, 7.7401463253, 6.6861119506, 6.740754889, 7.1061654625, 7.8072849019, 6.7065466324, 7.4918296506, 7.2361367928, 7.8459826909, 8.2448229784, 12.8411396258, 12.6936196106, 11.985032865, 12.6729528341, 13.378352665, 13.1446822034, 13.1254440465, 13.2037164414, 12.9402222152, 13.466498446, 12.9024963514, 13.8857598862, 13.2576750354, 9.2781931791, 8.9382132998, 9.0884695656, 9.080461612, 9.0903272779, 9.8286044069, 10.4947401005, 10.6659360514, 9.1369780474, 9.8052880973, 9.5229691492, 10.7971712442, 9.162506677, 12.4995506953, 12.5152493103, 12.7927852785, 13.0303983721, 12.109103588, 13.1907709613, 13.2980629749, 12.7312223381, 14.3319836381, 13.8856856732, 12.9429105173, 13.6234032462, 14.9871419763, 11.4167070957, 12.0622443156, 11.5288290421, 12.7950848055, 12.2018105064, 12.300439657, 12.3741974889, 12.2215291205, 13.2058303327, 13.4778966329, 12.8716265451, 12.4938346021, 13.9917171333, 17.6261181394, 18.0306168464, 18.4613922389, 19.1324268388, 19.1917842883, 18.7270624105, 19.0520425589, 18.4280677509, 19.459420989, 18.7670794726, 18.8027245471, 19.2222287854, 18.9300055987, 11.5079131779, 11.0922487521, 11.6658567479, 10.858033816, 9.977366855, 12.880651091, 11.4479959218, 11.7528732296, 11.9295317438, 11.7026503736, 12.67898185, 11.6577854559, 11.3264234768, 13.7453700883, 14.8580708601, 13.9674362984, 14.5370920116, 14.7499800173, 14.0391991165, 14.560852081, 13.5000562869, 12.9911620457, 15.8540447024, 15.0560012133, 15.0023763083, 14.6526322112, 14.1552037443, 13.2495794367, 14.1266577346, 13.6631412804, 14.3095479369, 14.2696488631, 14.0248507859, 13.6451859344, 13.8540790367, 14.6017639896, 14.9257345988, 14.9288799209, 14.7800388018, 14.1624225519, 13.9704424775, 15.3304792663, 13.7287623769, 14.4093659852, 14.1579161453, 15.2119878012, 14.1736251862, 15.0723855661, 14.9996740159, 15.043764438, 16.2575935502, 15.6167561526, 10.2121454476, 10.7071628622, 11.2630032475, 10.1115433787, 10.4303691912, 11.450933045, 11.2357913452, 11.2836421946, 12.1603326216, 10.1561960411, 10.9837703547, 10.5524871865, 11.0496650534, 11.8019159442, 12.2668266019, 11.4559698437, 10.6889619541, 12.1574167835, 11.77035963, 12.2234819278, 13.1912496553, 11.6959041846, 12.0218821641, 11.5847917678, 12.1508212257, 12.3653691566, 14.0717834031, 13.9017194733, 13.6896818335, 12.8219381465, 13.9515712037, 13.4878901181, 13.5663366903, 14.7698866956, 13.7459506981, 14.2057001753, 14.9239745639, 15.0617567635, 15.3721118258, 8.0518854599, 8.1590017681, 8.278269688, 9.0193685781, 8.6624064498, 8.0299934392, 8.129057928, 8.3905914418, 8.4535556206, 8.2061123566, 8.9326235245, 9.7314166211, 10.0804246162, 16.2545664928, 17.2572109812, 17.4349414991, 17.6863866397, 16.894621256, 16.9744888549, 18.0702391393, 17.026968762, 17.9778157216, 17.0473679201, 18.2342288525, 17.6813282932, 17.7675643621, 15.3769516713, 14.0851398882, 14.6143565691, 14.7892799275, 15.1327581575, 14.7696599461, 15.2400045477, 15.2241965765, 15.5581598562, 15.8517787282, 15.5225873145, 16.3280758601, 16.0059029786, 13.5135665994, 11.7238938271, 13.0432181346, 13.472584571, 13.7185317486, 12.59047135, 14.0476208463, 13.6663054684, 12.7798481751, 13.978531447, 14.2300761276, 13.0276658509, 14.7660788106, 12.2899358223, 12.4941017173, 11.5258893465, 11.7561649108, 12.9409325034, 11.4894195968, 12.5272849402, 13.3026607982, 12.176517092, 12.4459536149, 13.0792254017, 13.3333969434, 12.3660994776, 12.8951280351, 11.9873940936, 12.3552711049, 12.6162664943, 11.6312331335, 12.621878607, 11.7706082416, 12.6151521192, 12.8868988012, 12.4956213065, 13.0783868081, 12.9740483638, 12.8536843573, 10.5322884982, 10.3095656168, 10.5958178202, 9.5446471933, 9.5116621718, 10.7478545292, 10.9106163212, 9.9430070016, 11.3666842845, 10.6341994205, 11.1796926464, 11.8299400287, 11.2599998495, 10.3331338907, 10.8786246056, 10.361961281, 10.3180008811, 11.1949752298, 10.660986253, 11.4851093893, 10.3382989987, 11.7282986672, 11.9182559595, 11.5329533877, 11.7648480374, 12.0467202617, 13.19799333, 12.7196476155, 13.882500464, 13.6153658669, 13.4595949409, 13.5451076472, 13.8112491665, 14.6457010289, 14.2845283257, 14.3562295142, 15.5632444449, 14.1560946111, 14.4293750419] + }, + "params": { + "pattern": "single_baseline_multi_path_non_binary", + "n_switcher_groups": 78, + "n_realized_groups": 118, + "n_periods": 13, + "seed": 119, + "effects": 3, + "placebo": 1, + "by_path": 3, + "ci_level": 95 + }, + "results": { + "by_path": [ + { + "path": "0,1,1,1", + "frequency_rank": 1, + "horizons": { + "1": { + "effect": 1.4178257704, + "se": 0.15332691296, + "ci_lo": 1.1173105432, + "ci_hi": 1.7183409977, + "n_switchers": 32, + "n_obs": 178 + }, + "2": { + "effect": 1.1888829872, + "se": 0.13389814869, + "ci_lo": 0.9264474382, + "ci_hi": 1.4513185363, + "n_switchers": 32, + "n_obs": 150 + }, + "3": { + "effect": 1.1343468218, + "se": 0.13441857307, + "ci_lo": 0.87089125968, + "ci_hi": 1.3978023838, + "n_switchers": 32, + "n_obs": 124 + }, + "-1": { + "effect": -0.07428441879, + "se": 0.13418342659, + "ci_lo": -0.33727910222, + "ci_hi": 0.18871026464, + "n_switchers": 32, + "n_obs": 178 + } + } + }, + { + "path": "0,2,2,2", + "frequency_rank": 2, + "horizons": { + "1": { + "effect": 2.7581718804, + "se": 0.11089003798, + "ci_lo": 2.5408313997, + "ci_hi": 2.9755123611, + "n_switchers": 26, + "n_obs": 118 + }, + "2": { + "effect": 2.7158629099, + "se": 0.13743908945, + "ci_lo": 2.4464872445, + "ci_hi": 2.9852385753, + "n_switchers": 26, + "n_obs": 94 + }, + "3": { + "effect": 2.5576072535, + "se": 0.16260901373, + "ci_lo": 2.2388994431, + "ci_hi": 2.876315064, + "n_switchers": 26, + "n_obs": 74 + }, + "-1": { + "effect": -0.0033722906152, + "se": 0.13786455862, + "ci_lo": -0.27358186026, + "ci_hi": 0.26683727903, + "n_switchers": 26, + "n_obs": 118 + } + } + }, + { + "path": "0,1,2,2", + "frequency_rank": 3, + "horizons": { + "1": { + "effect": 1.2765396037, + "se": 0.19144937356, + "ci_lo": 0.90130572669, + "ci_hi": 1.6517734808, + "n_switchers": 20, + "n_obs": 68 + }, + "2": { + "effect": 2.9084357291, + "se": 0.20796144712, + "ci_lo": 2.5008387825, + "ci_hi": 3.3160326756, + "n_switchers": 20, + "n_obs": 60 + }, + "3": { + "effect": 2.8957739292, + "se": 0.19309948403, + "ci_lo": 2.517305895, + "ci_hi": 3.2742419633, + "n_switchers": 20, + "n_obs": 60 + }, + "-1": { + "effect": -0.10992248773, + "se": 0.17275607563, + "ci_lo": -0.44851817408, + "ci_hi": 0.22867319861, + "n_switchers": 20, + "n_obs": 68 + } + } + } + ] + } } }, "generator": "generate_reversible_did_data v1", diff --git a/diff_diff/chaisemartin_dhaultfoeuille.py b/diff_diff/chaisemartin_dhaultfoeuille.py index 2e6fe940..51c4339a 100644 --- a/diff_diff/chaisemartin_dhaultfoeuille.py +++ b/diff_diff/chaisemartin_dhaultfoeuille.py @@ -26,7 +26,7 @@ """ import warnings -from typing import Any, Dict, List, Optional, Tuple +from typing import Any, Dict, List, Optional, Sequence, Tuple import numpy as np import pandas as pd @@ -282,6 +282,52 @@ def _validate_and_aggregate_to_cells( return cell +def _validate_paths_of_interest( + paths_of_interest: Any, +) -> List[Tuple[int, ...]]: + """Validate and canonicalize ``paths_of_interest`` to ``List[Tuple[int, ...]]``. + + Rejects non-sequence inputs, empty lists, non-tuple/list path entries, + empty path entries, non-int elements (including ``bool`` and + ``np.bool_``), and entries with mixed lengths. Numpy integer types + (``np.integer``) are accepted and canonicalized to Python ``int`` + so the resulting tuples are usable as dict keys interchangeably with + paths emitted by ``_enumerate_treatment_paths`` (which casts via + ``int(round(float(v)))``). + """ + if not isinstance(paths_of_interest, (list, tuple)): + raise ValueError( + f"paths_of_interest must be a list/tuple of int tuples, " + f"got {type(paths_of_interest).__name__}." + ) + if len(paths_of_interest) == 0: + raise ValueError("paths_of_interest must be non-empty.") + canonical: List[Tuple[int, ...]] = [] + for i, p in enumerate(paths_of_interest): + if not isinstance(p, (list, tuple)): + raise ValueError( + f"paths_of_interest[{i}] must be a tuple/list of ints, " f"got {type(p).__name__}." + ) + if len(p) == 0: + raise ValueError(f"paths_of_interest[{i}] must be non-empty.") + canonical_path: List[int] = [] + for j, v in enumerate(p): + if isinstance(v, (bool, np.bool_)) or not isinstance(v, (int, np.integer)): + raise ValueError( + f"paths_of_interest[{i}][{j}] must be an int, got " + f"{v!r} of type {type(v).__name__}." + ) + canonical_path.append(int(v)) + canonical.append(tuple(canonical_path)) + lens = {len(p) for p in canonical} + if len(lens) > 1: + raise ValueError( + f"paths_of_interest entries must all have the same length " + f"(L_max+1); got mixed lengths {sorted(lens)}." + ) + return canonical + + # ============================================================================= # Main estimator class # ============================================================================= @@ -320,8 +366,10 @@ class ChaisemartinDHaultfoeuille(ChaisemartinDHaultfoeuilleBootstrapMixin): HonestDiD sensitivity integration on placebos via ``honest_did=True`` - Per-path event-study disaggregation via ``by_path=k`` (top-k most common observed treatment paths within the window - ``[F_g-1, F_g-1+L_max]``; requires ``drop_larger_lower=False`` and - binary treatment) + ``[F_g-1, F_g-1+L_max]``; requires ``drop_larger_lower=False``; + supports binary or integer-coded discrete treatment) or via + ``paths_of_interest=[(...), ...]`` for an explicit user-specified + path subset (Python-only API; mutex with ``by_path=k``) - Survey support via ``survey_design=``: pweight with strata/PSU/FPC via Taylor Series Linearization (analytical) or replicate-weight variance (BRR/Fay/JK1/JKn/SDR) @@ -406,12 +454,19 @@ class ChaisemartinDHaultfoeuille(ChaisemartinDHaultfoeuilleBootstrapMixin): Requires ``drop_larger_lower=False`` (multi-switch groups are the object of interest) and ``L_max >= 1`` (the path window - depends on ``L_max``). Binary treatment only — non-binary - treatment + ``by_path`` is deferred. Also incompatible with - ``heterogeneity``, ``design2``, ``honest_did``, and - ``survey_design`` (each combination raises + depends on ``L_max``). Compatible with non-binary integer-coded + treatment (D in Z); path tuples become integer-state tuples + like ``(0, 2, 2, 2)``. D values must be integer-valued + (``D == round(D)``); a ``ValueError`` is raised at fit-time on + continuous D. Incompatible with ``heterogeneity``, ``design2``, + ``honest_did``, and ``survey_design`` (each combination raises ``NotImplementedError`` in the current release). + Mutually exclusive with ``paths_of_interest`` — use + ``by_path=k`` for top-k automatic ranking by frequency, or + ``paths_of_interest=[(...), ...]`` for an explicit user- + specified path list. Setting both raises ``ValueError``. + Compatible with ``controls`` (DID^X residualization) -- the per-baseline OLS residualization runs once on first-differenced ``Y`` BEFORE path enumeration, so per-path point estimates, @@ -538,6 +593,53 @@ class ChaisemartinDHaultfoeuille(ChaisemartinDHaultfoeuilleBootstrapMixin): by the path tuple, with nested ``"horizons"`` dicts per horizon ``l``. Also available via ``results.to_dataframe(level="by_path")``. + paths_of_interest : list of tuple of int, optional, default=None + Explicit user-specified treatment paths to disaggregate by, as + an alternative to ``by_path=k``'s top-k automatic ranking. + Each path tuple must have length ``L_max + 1`` and represents + the treatment trajectory in the window + ``[F_g - 1, F_g, ..., F_g - 1 + L_max]``, e.g. + ``[(0, 1, 1, 1), (0, 1, 0, 0)]`` for two paths under + ``L_max=3``. Mutually exclusive with ``by_path``; setting both + raises ``ValueError``. + + Validation: + + - Each path element must be an ``int`` (``bool`` and + ``np.bool_`` rejected; ``np.integer`` accepted and + canonicalized to Python ``int``). + - All paths must have the same length (uniformity validated + at ``__init__``; length match against ``L_max + 1`` + validated at fit-time). + - Empty list raises ``ValueError``. + - Duplicate paths are deduplicated with a ``UserWarning``. + - A path with zero observed groups in the panel emits a + ``UserWarning`` and is omitted from ``path_effects``. + + Compatible with non-binary integer treatment (paths can + contain integer states like ``(0, 2, 2)``). + + Compatible with all downstream surfaces inherited by + ``by_path``: bootstrap, per-path placebos, per-path joint + sup-t bands, ``controls``, ``trends_linear``, + ``trends_nonparam``. Mechanical extension to path + enumeration; no methodology change. + + **Order semantics**: paths appear in + ``results.path_effects`` in the user-specified order, modulo + deduplication and unobserved-path filtering. + + **Python-only API extension; no R equivalent.** R's + ``did_multiplegt_dyn(..., by_path=k)`` only accepts a positive + int (top-k) or ``-1`` (all paths); there is no list-based + path selection in R. + + Results expose the same surfaces as ``by_path``: + ``results.path_effects`` (dict keyed by path tuple), + ``results.path_placebo_event_study``, + ``results.path_sup_t_bands``, + ``results.path_cumulated_event_study`` (under + ``trends_linear``), and the ``level="by_path"`` DataFrame. rank_deficient_action : str, default="warn" Action when the TWFE decomposition diagnostic OLS encounters a rank-deficient design matrix: ``"warn"``, ``"error"``, or @@ -584,6 +686,7 @@ def __init__( twfe_diagnostic: bool = True, drop_larger_lower: bool = True, by_path: Optional[int] = None, + paths_of_interest: Optional[Sequence[Sequence[int]]] = None, rank_deficient_action: str = "warn", ) -> None: # Parameter validation @@ -610,9 +713,18 @@ def __init__( if by_path <= 0: raise ValueError( f"by_path must be a positive int (top-k most common paths), " - f"got {by_path}. Use by_path=None to disable; explicit path " - f"selection via paths_of_interest is a future extension." + f"got {by_path}. Use by_path=None to disable, or " + f"paths_of_interest for explicit path selection." ) + if paths_of_interest is not None: + paths_of_interest = _validate_paths_of_interest(paths_of_interest) + if by_path is not None and paths_of_interest is not None: + raise ValueError( + "by_path and paths_of_interest are mutually exclusive. " + "Use by_path=k for top-k automatic ranking, OR " + "paths_of_interest=[(...), ...] for explicit user-" + "specified paths. Set one and leave the other as None." + ) if cluster is not None: raise NotImplementedError( f"cluster={cluster!r}: user-specified clustering is not " @@ -637,6 +749,7 @@ def __init__( self.twfe_diagnostic = twfe_diagnostic self.drop_larger_lower = drop_larger_lower self.by_path = by_path + self.paths_of_interest = paths_of_interest self.rank_deficient_action = rank_deficient_action self.is_fitted_ = False @@ -658,6 +771,7 @@ def get_params(self) -> Dict[str, Any]: "twfe_diagnostic": self.twfe_diagnostic, "drop_larger_lower": self.drop_larger_lower, "by_path": self.by_path, + "paths_of_interest": self.paths_of_interest, "rank_deficient_action": self.rank_deficient_action, } @@ -665,14 +779,30 @@ def set_params(self, **params: Any) -> "ChaisemartinDHaultfoeuille": """ Set estimator parameters (sklearn-compatible). - Re-runs the same validation rules as ``__init__`` so invalid - parameter combinations cannot be introduced after construction. + **Transactional**: validation runs after the candidate mutations, + and if any rule fails the estimator state is rolled back to its + pre-call values before the exception is re-raised. Callers can + therefore retry with corrected params on the same instance + without repairing inconsistent intermediate state. """ - for key, value in params.items(): + # Snapshot current values for the keys we are about to set so + # we can roll back on validation failure (transactional semantics). + for key in params: if not hasattr(self, key): raise ValueError(f"Unknown parameter: {key}") - setattr(self, key, value) + snapshot = {key: getattr(self, key) for key in params} + try: + for key, value in params.items(): + setattr(self, key, value) + self._validate_invariants() + except Exception: + for key, value in snapshot.items(): + setattr(self, key, value) + raise + return self + def _validate_invariants(self) -> None: + """Run the post-mutation validation rules. Mirrors `__init__`.""" # Re-run __init__ validation rules so the post-set state is valid. if self.rank_deficient_action not in ("warn", "error", "silent"): raise ValueError( @@ -697,9 +827,18 @@ def set_params(self, **params: Any) -> "ChaisemartinDHaultfoeuille": if self.by_path <= 0: raise ValueError( f"by_path must be a positive int (top-k most common paths), " - f"got {self.by_path}. Use by_path=None to disable; explicit " - f"path selection via paths_of_interest is a future extension." + f"got {self.by_path}. Use by_path=None to disable, or " + f"paths_of_interest for explicit path selection." ) + if self.paths_of_interest is not None: + self.paths_of_interest = _validate_paths_of_interest(self.paths_of_interest) + if self.by_path is not None and self.paths_of_interest is not None: + raise ValueError( + "by_path and paths_of_interest are mutually exclusive. " + "Use by_path=k for top-k automatic ranking, OR " + "paths_of_interest=[(...), ...] for explicit user-" + "specified paths. Set one and leave the other as None." + ) if self.cluster is not None: raise NotImplementedError( f"cluster={self.cluster!r}: user-specified clustering is " @@ -711,7 +850,6 @@ def set_params(self, **params: Any) -> "ChaisemartinDHaultfoeuille": f"clustering is reserved for a future phase. See REGISTRY.md " f"ChaisemartinDHaultfoeuille section for the full contract." ) - return self # ------------------------------------------------------------------ # fit @@ -1041,44 +1179,63 @@ def fit( ) # ------------------------------------------------------------------ - # by_path preconditions and Phase 3 compatibility gates + # by_path / paths_of_interest preconditions and Phase 3 + # compatibility gates # ------------------------------------------------------------------ - if self.by_path is not None: + if self.by_path is not None or self.paths_of_interest is not None: if self.drop_larger_lower: raise ValueError( - "by_path requires drop_larger_lower=False because " - "multi-switch groups are the object of interest for " - "per-path disaggregation, but the default " - "drop_larger_lower=True filter removes them. Construct " - "the estimator with " + "by_path / paths_of_interest requires " + "drop_larger_lower=False because multi-switch groups " + "are the object of interest for per-path " + "disaggregation, but the default " + "drop_larger_lower=True filter removes them. " + "Construct the estimator with " "ChaisemartinDHaultfoeuille(drop_larger_lower=False, " - "by_path=k)." + "by_path=k) or " + "ChaisemartinDHaultfoeuille(" + "drop_larger_lower=False, " + "paths_of_interest=[(...), ...])." ) if L_max is None: raise ValueError( - "by_path requires L_max >= 1. The path window spans " - "[F_g - 1, F_g - 1 + L_max] and therefore depends on " - "the event-study horizon. Set L_max when calling fit()." + "by_path / paths_of_interest requires L_max >= 1. " + "The path window spans [F_g - 1, F_g - 1 + L_max] " + "and therefore depends on the event-study horizon. " + "Set L_max when calling fit()." ) + if self.paths_of_interest is not None: + expected_len = L_max + 1 + for p in self.paths_of_interest: + if len(p) != expected_len: + raise ValueError( + f"paths_of_interest entries must have " + f"length L_max+1={expected_len} (window " + f"[F_g-1, ..., F_g-1+L_max]); got path " + f"{p!r} of length {len(p)}." + ) if heterogeneity is not None: raise NotImplementedError( - "by_path combined with heterogeneity testing is " - "deferred to a future release." + "by_path / paths_of_interest combined with " + "heterogeneity testing is deferred to a future release." ) if design2: raise NotImplementedError( - "by_path combined with design2 is deferred to a future " "release." + "by_path / paths_of_interest combined with design2 " + "is deferred to a future release." ) if honest_did: raise NotImplementedError( - "by_path combined with honest_did (HonestDiD sensitivity " - "analysis) is deferred to a future release." + "by_path / paths_of_interest combined with honest_did " + "(HonestDiD sensitivity analysis) is deferred to a " + "future release." ) if survey_design is not None: raise NotImplementedError( - "by_path combined with survey_design is deferred to a " - "future release: the cell-period IF allocator under " - "path subsets has not been derived." + "by_path / paths_of_interest combined with " + "survey_design is deferred to a future release: the " + "cell-period IF allocator under path subsets has not " + "been derived." ) # ------------------------------------------------------------------ @@ -1574,18 +1731,19 @@ def fit( # `D_{g,1}` values across the never-treated / always-treated # control mix). SE inheritance (cross-path cohort-sharing) is # documented separately in REGISTRY.md. - if self.by_path is not None: + if self.by_path is not None or self.paths_of_interest is not None: _switcher_mask = first_switch_idx_arr >= 0 if _switcher_mask.any(): _switcher_baselines = baselines[_switcher_mask] if np.unique(_switcher_baselines).size > 1: warnings.warn( - "by_path + controls: switcher baselines D_{g,1} " - "take multiple values in this panel. Python " - "residualizes once on the full panel before path " - "enumeration; R `did_multiplegt_dyn(..., by_path, " - "controls)` re-runs residualization per path on " - "the path-restricted subsample, so per-path point " + "by_path / paths_of_interest + controls: " + "switcher baselines D_{g,1} take multiple values " + "in this panel. Python residualizes once on the " + "full panel before path enumeration; R " + "`did_multiplegt_dyn(..., by_path, controls)` " + "re-runs residualization per path on the " + "path-restricted subsample, so per-path point " "estimates can diverge between Python and R on " "this panel. See `docs/methodology/REGISTRY.md` " "(`Note (Phase 3 by_path ...)` -> Per-path " @@ -1678,16 +1836,16 @@ def fit( # that disagree with R. The check filters to switcher groups # only (never-switchers / always-treated controls don't # contribute to switcher baseline multiplicity). - if self.by_path is not None: + if self.by_path is not None or self.paths_of_interest is not None: _switcher_mask_tl = first_switch_idx_arr >= 0 if _switcher_mask_tl.any(): _switcher_baselines_tl = baselines[_switcher_mask_tl] if np.unique(_switcher_baselines_tl).size > 1: warnings.warn( - "by_path + trends_linear: switcher baselines " - "D_{g,1} take multiple values in this panel. " - "Python first-differences once on the full " - "panel before path enumeration; R " + "by_path / paths_of_interest + trends_linear: " + "switcher baselines D_{g,1} take multiple values " + "in this panel. Python first-differences once on " + "the full panel before path enumeration; R " "`did_multiplegt_dyn(..., by_path, trends_lin)` " "re-runs the full pipeline (including " "first-differencing) on each path's restricted " @@ -1719,8 +1877,9 @@ def fit( _f_g_three_count = int((first_switch_idx_arr == 2).sum()) if _f_g_three_count > 0: warnings.warn( - f"by_path + trends_linear: {_f_g_three_count} " - f"switching group(s) have F_g=3 (exactly 2 " + f"by_path / paths_of_interest + trends_linear: " + f"{_f_g_three_count} switching group(s) have " + f"F_g=3 (exactly 2 " f"pre-switch periods). After first-differencing " f"and the time==1 filter, these groups have " f"only 1 valid pre-window Z value, which " @@ -1867,13 +2026,18 @@ def fit( "use the per-group DID_{g,l} building block which handles " "non-binary treatment." ) - if self.by_path is not None and not is_binary: - raise NotImplementedError( - "by_path combined with non-binary treatment is deferred to " - "a future release. Path enumeration requires integer-valued " - "treatment states to construct deterministic path tuples. " - "Use by_path=None with non-binary treatment." - ) + if (self.by_path is not None or self.paths_of_interest is not None) and not is_binary: + finite_D = D_mat[~np.isnan(D_mat)] + if finite_D.size > 0 and not np.all(finite_D == np.round(finite_D)): + bad_examples = np.unique(finite_D[finite_D != np.round(finite_D)])[:3] + raise ValueError( + f"by_path / paths_of_interest with non-binary " + f"treatment requires integer-coded treatment values " + f"(D in Z). Found non-integer values: " + f"{bad_examples.tolist()!r}. Round/discretize D " + f"before fitting, or set by_path=None and " + f"paths_of_interest=None with continuous treatment." + ) if N_S == 0 and (L_max is None or is_binary): raise ValueError( "No switching cells found in the data after filtering: every " @@ -2160,11 +2324,11 @@ def fit( # by_path disaggregation by observed treatment trajectory path_effects: Optional[Dict[Tuple[int, ...], Dict[str, Any]]] = None path_placebos: Optional[Dict[Tuple[int, ...], Dict[int, Dict[str, Any]]]] = None - path_cumulated_event_study: Optional[ - Dict[Tuple[int, ...], Dict[int, Dict[str, Any]]] - ] = None + path_cumulated_event_study: Optional[Dict[Tuple[int, ...], Dict[int, Dict[str, Any]]]] = ( + None + ) if ( - self.by_path is not None + (self.by_path is not None or self.paths_of_interest is not None) and L_max is not None and L_max >= 1 and multi_horizon_dids is not None @@ -2180,6 +2344,7 @@ def fit( T_g=T_g_arr, L_max=L_max, by_path=self.by_path, + paths_of_interest=self.paths_of_interest, eligible_mask_var=eligible_mask_var, multi_horizon_dids=multi_horizon_dids, all_groups=all_groups, @@ -2330,7 +2495,11 @@ def fit( # path cohort-sharing SE deviation from R documented for # `path_effects` (full-panel cohort-centered plug-in vs # R's per-path re-run). - if self.by_path is not None and self.placebo and multi_horizon_placebos is not None: + if ( + (self.by_path is not None or self.paths_of_interest is not None) + and self.placebo + and multi_horizon_placebos is not None + ): _df_s_bp_pl = _effective_df_survey(resolved_survey, _replicate_n_valid_list) path_placebos = _compute_path_placebos( D_mat=D_mat, @@ -2342,6 +2511,7 @@ def fit( T_g=T_g_arr, L_max=L_max, by_path=self.by_path, + paths_of_interest=self.paths_of_interest, eligible_mask_var=eligible_mask_var, multi_horizon_placebos=multi_horizon_placebos, alpha=self.alpha, @@ -2870,7 +3040,7 @@ def fit( # architectural preference — see `_collect_path_bootstrap_inputs`). path_bootstrap_inputs = None if ( - self.by_path is not None + (self.by_path is not None or self.paths_of_interest is not None) and L_max is not None and L_max >= 1 and multi_horizon_dids is not None @@ -2887,6 +3057,7 @@ def fit( T_g=T_g_arr, L_max=L_max, by_path=self.by_path, + paths_of_interest=self.paths_of_interest, eligible_mask_var=eligible_mask_var, multi_horizon_dids=multi_horizon_dids, path_effects=path_effects, @@ -2900,7 +3071,7 @@ def fit( # empty dict. path_placebo_bootstrap_inputs = None if ( - self.by_path is not None + (self.by_path is not None or self.paths_of_interest is not None) and self.placebo and L_max is not None and L_max >= 1 @@ -2918,6 +3089,7 @@ def fit( T_g=T_g_arr, L_max=L_max, by_path=self.by_path, + paths_of_interest=self.paths_of_interest, eligible_mask_var=eligible_mask_var, multi_horizon_placebos=multi_horizon_placebos, path_placebos=path_placebos, @@ -3177,7 +3349,7 @@ def fit( # cumulated SE / t / p / CI, regardless of whether the source # was an analytical singularity or a non-finite bootstrap draw. if ( - self.by_path is not None + (self.by_path is not None or self.paths_of_interest is not None) and _is_trends_linear and L_max is not None and L_max >= 1 @@ -3185,9 +3357,7 @@ def fit( and path_effects is not None and len(path_effects) > 0 ): - _df_s_bp_cum = _effective_df_survey( - resolved_survey, _replicate_n_valid_list - ) + _df_s_bp_cum = _effective_df_survey(resolved_survey, _replicate_n_valid_list) path_cumulated_event_study = _compute_path_cumulated_event_study( D_mat=D_mat, N_mat=N_mat, @@ -3195,6 +3365,7 @@ def fit( switch_direction=switch_direction_arr, L_max=L_max, by_path=self.by_path, + paths_of_interest=self.paths_of_interest, multi_horizon_dids=multi_horizon_dids, path_effects=path_effects, alpha=self.alpha, @@ -3945,7 +4116,10 @@ def fit( ).items() if np.isfinite(crit) } - if (self.by_path is not None and self.n_bootstrap > 0) + if ( + (self.by_path is not None or self.paths_of_interest is not None) + and self.n_bootstrap > 0 + ) else None ), survey_metadata=survey_metadata, @@ -5459,14 +5633,16 @@ def _enumerate_treatment_paths( first_switch_idx: np.ndarray, N_mat: np.ndarray, L_max: int, - by_path: int, + by_path: Optional[int], + paths_of_interest: Optional[List[Tuple[int, ...]]] = None, ) -> Tuple[ List[Tuple[int, ...]], Dict[Tuple[int, ...], np.ndarray], Dict[Tuple[int, ...], int], ]: """ - Enumerate observed treatment paths and select the top-``by_path`` most common. + Enumerate observed treatment paths and select either the top-``by_path`` + most common or the user-specified ``paths_of_interest`` subset. For each switcher group ``g``, the path is the treatment tuple ``(D_{g, F_g-1}, D_{g, F_g}, ..., D_{g, F_g-1+L_max})`` — length @@ -5479,24 +5655,37 @@ def _enumerate_treatment_paths( the number of observed paths, all observed paths are returned with a ``UserWarning``. + When ``paths_of_interest`` is provided, the user-specified subset is + used instead of the top-k ranking. Duplicate paths emit a + ``UserWarning`` and are deduplicated; paths not observed in the + panel emit a ``UserWarning`` and are omitted from the result. + Parameters ---------- D_mat : np.ndarray of shape (n_groups, n_periods) - Treatment matrix. Must be binary (0/1); non-binary treatment is - gated out upstream. + Treatment matrix. Binary (0/1) or integer-coded discrete + treatment (D in Z); upstream validation enforces D == round(D) + when ``not is_binary``. first_switch_idx : np.ndarray of shape (n_groups,) Index of first switch per group; ``-1`` for never-switching groups. N_mat : np.ndarray of shape (n_groups, n_periods) Cell-count matrix (zero where the cell is unobserved). L_max : int Event-study horizon; window length is ``L_max + 1``. - by_path : int - Number of most-common paths to select. + by_path : int or None + Number of most-common paths to select. Mutually exclusive with + ``paths_of_interest``; exactly one of the two is non-None when + the per-path branch fires. + paths_of_interest : list of tuple[int, ...] or None, default None + User-specified path subset. When provided, ``by_path`` is + ignored. Returns ------- selected_paths : list of tuple[int, ...] - Selected path tuples, ordered by descending frequency. + Selected path tuples. Under ``by_path``, ordered by descending + frequency. Under ``paths_of_interest``, in user-specified order + modulo deduplication and unobserved-path filtering. path_to_group_mask : dict[tuple[int, ...], np.ndarray of bool] Per-path boolean mask over all ``n_groups`` identifying switchers that follow that path. @@ -5525,20 +5714,50 @@ def _enumerate_treatment_paths( for path in path_of_group.values(): path_counts[path] = path_counts.get(path, 0) + 1 - observed_paths = sorted(path_counts.keys(), key=lambda p: (-path_counts[p], p)) - n_observed = len(observed_paths) - - if by_path >= n_observed: - if by_path > n_observed and n_observed > 0: - warnings.warn( - f"by_path={by_path} exceeds the number of observed paths " - f"({n_observed}). Returning all observed paths.", - UserWarning, - stacklevel=2, - ) - selected_paths = observed_paths + if paths_of_interest is not None: + # Canonicalization (Tuple[int, ...] with Python int) happens in + # __init__ / set_params via _validate_paths_of_interest, so + # duplicates such as `[(np.int64(0), 1, 1, 1), (0, 1, 1, 1)]` + # collapse to the same tuple here and the seen-set check fires. + observed_paths_set = set(path_counts.keys()) + seen: set = set() + selected_paths: List[Tuple[int, ...]] = [] + for p in paths_of_interest: + if p in seen: + warnings.warn( + f"paths_of_interest contains duplicate path {p!r}; " f"deduplicating.", + UserWarning, + stacklevel=2, + ) + continue + seen.add(p) + if p not in observed_paths_set: + warnings.warn( + f"paths_of_interest path {p!r} has zero observed " + f"groups in the panel; this path will be omitted " + f"from path_effects.", + UserWarning, + stacklevel=2, + ) + continue + selected_paths.append(p) else: - selected_paths = observed_paths[:by_path] + observed_paths = sorted(path_counts.keys(), key=lambda p: (-path_counts[p], p)) + n_observed = len(observed_paths) + if by_path is None: + # Defensive: caller always passes one of the two selectors. + selected_paths = [] + elif by_path >= n_observed: + if by_path > n_observed and n_observed > 0: + warnings.warn( + f"by_path={by_path} exceeds the number of observed " + f"paths ({n_observed}). Returning all observed paths.", + UserWarning, + stacklevel=2, + ) + selected_paths = observed_paths + else: + selected_paths = observed_paths[:by_path] path_to_group_mask: Dict[Tuple[int, ...], np.ndarray] = {} for path in selected_paths: @@ -5561,13 +5780,14 @@ def _compute_path_effects( switch_direction: np.ndarray, T_g: np.ndarray, L_max: int, - by_path: int, + by_path: Optional[int], eligible_mask_var: np.ndarray, multi_horizon_dids: Dict[int, Dict[str, Any]], all_groups: List[Any], alpha: float, df_inference: Optional[int] = None, set_ids: Optional[np.ndarray] = None, + paths_of_interest: Optional[List[Tuple[int, ...]]] = None, ) -> Optional[Dict[Tuple[int, ...], Dict[str, Any]]]: """ Compute per-path event-study effects using the joiners/leavers IF pattern. @@ -5600,19 +5820,38 @@ def _compute_path_effects( N_mat=N_mat, L_max=L_max, by_path=by_path, + paths_of_interest=paths_of_interest, ) if not selected_paths: - warnings.warn( - f"by_path={by_path} was requested but no observed treatment " - f"path has a complete window [F_g-1, F_g-1+L_max={L_max}] " - f"within the panel. results.path_effects is populated as an " - f"empty dict to signal 'requested but empty'. Extend the " - f"panel so switchers have L_max+1 consecutive observed cells " - f"starting at F_g-1, or reduce L_max.", - UserWarning, - stacklevel=2, - ) + if paths_of_interest is not None: + # Every requested path was unobserved (each emitted its own + # per-path "zero observed groups" warning inside the + # enumerator). Distinguish from the by_path=k case where + # the panel itself has no complete-window path. + warnings.warn( + "paths_of_interest was requested but every " + "user-specified path was either unobserved in the " + "panel or had a window outside the L_max+1 " + "convention (per-path 'zero observed groups' " + "UserWarnings already issued). results.path_effects " + "is populated as an empty dict to signal 'requested " + "but empty'.", + UserWarning, + stacklevel=2, + ) + else: + warnings.warn( + f"by_path={by_path} was requested but no observed " + f"treatment path has a complete window [F_g-1, " + f"F_g-1+L_max={L_max}] within the panel. " + f"results.path_effects is populated as an empty dict " + f"to signal 'requested but empty'. Extend the panel " + f"so switchers have L_max+1 consecutive observed " + f"cells starting at F_g-1, or reduce L_max.", + UserWarning, + stacklevel=2, + ) return {} # Cohort ids for the variance-eligible set (same construction as the @@ -5640,7 +5879,20 @@ def _compute_path_effects( path_effects: Dict[Tuple[int, ...], Dict[str, Any]] = {} - for rank, path in enumerate(selected_paths, start=1): + # `frequency_rank` is the within-selected-paths rank by descending + # group count (lex tiebreak on the path tuple). Decoupled from the + # iteration order over `selected_paths` so that under + # `paths_of_interest` (user-specified order) the rank still + # reflects true frequency. Under `by_path=k`, `selected_paths` is + # already sorted by descending frequency so ranks coincide with + # iteration order. + rank_sorted_paths = sorted( + selected_paths, + key=lambda p: (-path_to_count[p], p), + ) + path_to_freq_rank = {p: i + 1 for i, p in enumerate(rank_sorted_paths)} + + for path in selected_paths: switcher_mask = path_to_group_mask[path] n_path_groups = int(switcher_mask.sum()) @@ -5729,7 +5981,7 @@ def _compute_path_effects( path_effects[path] = { "n_groups": n_path_groups, - "frequency_rank": rank, + "frequency_rank": path_to_freq_rank[path], "horizons": horizons, } @@ -5742,11 +5994,12 @@ def _compute_path_cumulated_event_study( first_switch_idx: np.ndarray, switch_direction: np.ndarray, L_max: int, - by_path: int, + by_path: Optional[int], multi_horizon_dids: Dict[int, Dict[str, Any]], path_effects: Dict[Tuple[int, ...], Dict[str, Any]], alpha: float, df_inference: Optional[int] = None, + paths_of_interest: Optional[List[Tuple[int, ...]]] = None, ) -> Dict[Tuple[int, ...], Dict[int, Dict[str, Any]]]: """ Per-path cumulated level effects under ``trends_linear=True``. @@ -5779,6 +6032,7 @@ def _compute_path_cumulated_event_study( N_mat=N_mat, L_max=L_max, by_path=by_path, + paths_of_interest=paths_of_interest, ) n_groups_total = D_mat.shape[0] @@ -5829,15 +6083,13 @@ def _compute_path_cumulated_event_study( } continue cum_effect = float( - np.sum(S_arr[eligible_path] * running_per_group[eligible_path]) - / n_l_path + np.sum(S_arr[eligible_path] * running_per_group[eligible_path]) / n_l_path ) # Conservative SE upper bound: sum of per-horizon component # SEs from path_effects (matches global formula at :3402-3413). # NaN-consistency: any non-finite component SE -> cumulated NaN. component_ses = [ - path_horizons_anal.get(ll, {}).get("se", float("nan")) - for ll in range(1, l_h + 1) + path_horizons_anal.get(ll, {}).get("se", float("nan")) for ll in range(1, l_h + 1) ] if all(np.isfinite(s) for s in component_ses): cum_se = float(sum(component_ses)) @@ -5872,12 +6124,13 @@ def _compute_path_placebos( switch_direction: np.ndarray, T_g: np.ndarray, L_max: int, - by_path: int, + by_path: Optional[int], eligible_mask_var: np.ndarray, multi_horizon_placebos: Dict[int, Dict[str, Any]], alpha: float, df_inference: Optional[int] = None, set_ids: Optional[np.ndarray] = None, + paths_of_interest: Optional[List[Tuple[int, ...]]] = None, ) -> Optional[Dict[Tuple[int, ...], Dict[int, Dict[str, Any]]]]: """ Compute per-path backward-horizon placebos ``DID^{pl}_{path, l}``. @@ -5926,6 +6179,7 @@ def _compute_path_placebos( N_mat=N_mat, L_max=L_max, by_path=by_path, + paths_of_interest=paths_of_interest, ) if not selected_paths: @@ -6046,11 +6300,12 @@ def _collect_path_bootstrap_inputs( switch_direction: np.ndarray, T_g: np.ndarray, L_max: int, - by_path: int, + by_path: Optional[int], eligible_mask_var: np.ndarray, multi_horizon_dids: Dict[int, Dict[str, Any]], path_effects: Dict[Tuple[int, ...], Dict[str, Any]], set_ids: Optional[np.ndarray] = None, + paths_of_interest: Optional[List[Tuple[int, ...]]] = None, ) -> Dict[Tuple[int, ...], Dict[int, Tuple[np.ndarray, int, float, None]]]: """ Collect per-(path, horizon) inputs for the bootstrap mixin. @@ -6092,6 +6347,7 @@ def _collect_path_bootstrap_inputs( N_mat=N_mat, L_max=L_max, by_path=by_path, + paths_of_interest=paths_of_interest, ) n_groups = D_mat.shape[0] @@ -6171,11 +6427,12 @@ def _collect_path_placebo_bootstrap_inputs( switch_direction: np.ndarray, T_g: np.ndarray, L_max: int, - by_path: int, + by_path: Optional[int], eligible_mask_var: np.ndarray, multi_horizon_placebos: Dict[int, Dict[str, Any]], path_placebos: Dict[Tuple[int, ...], Dict[int, Dict[str, Any]]], set_ids: Optional[np.ndarray] = None, + paths_of_interest: Optional[List[Tuple[int, ...]]] = None, ) -> Dict[Tuple[int, ...], Dict[int, Tuple[np.ndarray, int, float, None]]]: """ Collect per-(path, lag) inputs for the placebo bootstrap mixin @@ -6215,6 +6472,7 @@ def _collect_path_placebo_bootstrap_inputs( N_mat=N_mat, L_max=L_max, by_path=by_path, + paths_of_interest=paths_of_interest, ) n_groups = D_mat.shape[0] @@ -7928,17 +8186,12 @@ def chaisemartin_dhaultfoeuille( ------- ChaisemartinDHaultfoeuilleResults """ + import inspect + init_keys = { - "alpha", - "cluster", - "n_bootstrap", - "bootstrap_weights", - "seed", - "placebo", - "twfe_diagnostic", - "drop_larger_lower", - "by_path", - "rank_deficient_action", + name + for name, p in inspect.signature(ChaisemartinDHaultfoeuille.__init__).parameters.items() + if p.kind not in (p.VAR_POSITIONAL, p.VAR_KEYWORD) and name != "self" } init_kwargs = {k: v for k, v in kwargs.items() if k in init_keys} fit_kwargs = {k: v for k, v in kwargs.items() if k not in init_keys} diff --git a/diff_diff/chaisemartin_dhaultfoeuille_results.py b/diff_diff/chaisemartin_dhaultfoeuille_results.py index 6ddfb0d6..5fcd2e7d 100644 --- a/diff_diff/chaisemartin_dhaultfoeuille_results.py +++ b/diff_diff/chaisemartin_dhaultfoeuille_results.py @@ -394,10 +394,15 @@ class ChaisemartinDHaultfoeuilleResults: path_effects : dict, optional Per-path event-study effects keyed by observed treatment trajectory (tuple of int). Populated when ``by_path`` is a - positive int at estimator construction. Each entry holds + positive int OR ``paths_of_interest`` is a list of int tuples + at estimator construction. Each entry holds ``{"n_groups": int, "frequency_rank": int, "horizons": {l: {"effect", "se", "t_stat", "p_value", - "conf_int", "n_obs"}}}`` for ``l = 1..L_max``. + "conf_int", "n_obs"}}}`` for ``l = 1..L_max``. Under + ``paths_of_interest``, dict-insertion order matches the user- + specified path order; ``frequency_rank`` is the within- + selected-paths rank by descending observed-group count + (decoupled from iteration order). path_placebo_event_study : dict, optional Per-path backward-horizon placebos ``DID^{pl}_{path, l}`` for ``l = 1..L_max``, keyed by observed treatment trajectory (tuple @@ -407,11 +412,12 @@ class ChaisemartinDHaultfoeuilleResults: **path_placebo_event_study[p]}`` view is well-formed across forward and backward horizons. Each inner entry holds ``{"effect", "se", "t_stat", "p_value", "conf_int", "n_obs"}``. - Populated when ``by_path`` is a positive int AND - ``placebo=True`` AND ``L_max >= 1``. Empty-state contract - mirrors ``path_effects``: ``None`` when ``by_path + placebo`` - was not requested; ``{}`` when requested but no observed path - has a complete window ``[F_g-1, F_g-1+L_max]`` within the + Populated when (``by_path`` is a positive int OR + ``paths_of_interest`` is set) AND ``placebo=True`` AND + ``L_max >= 1``. Empty-state contract mirrors ``path_effects``: + ``None`` when ``by_path / paths_of_interest + placebo`` was + not requested; ``{}`` when requested but no observed path has + a complete window ``[F_g-1, F_g-1+L_max]`` within the panel (the same regime where ``path_effects`` returns ``{}``, with the same ``UserWarning`` at fit-time). Downstream callers should distinguish the two states. Inherits the cross-path @@ -424,9 +430,9 @@ class ChaisemartinDHaultfoeuilleResults: keyed by observed treatment trajectory (tuple of int). Inner dict is keyed by horizon directly (no ``"horizons"`` wrapper); each entry holds ``{"effect", "se", "t_stat", "p_value", - "conf_int", "n_obs"}``. Populated when ``by_path`` is a - positive int AND ``trends_linear=True`` AND ``L_max >= 1``; - ``None`` otherwise. Mirrors the global ``linear_trends_effects`` + "conf_int", "n_obs"}``. Populated when (``by_path`` is a + positive int OR ``paths_of_interest`` is set) AND + ``trends_linear=True`` AND ``L_max >= 1``; ``None`` otherwise. Mirrors the global ``linear_trends_effects`` cumulation: SE on the cumulated layer is the conservative upper bound (sum of per-horizon component SEs from ``path_effects[path]["horizons"][l]["se"]``, NaN-consistent). @@ -443,13 +449,15 @@ class ChaisemartinDHaultfoeuilleResults: observed treatment trajectory (tuple of int). Each entry holds ``{"crit_value": float, "alpha": float, "n_bootstrap": int, "method": str, "n_valid_horizons": int}``. Populated when - ``by_path`` is a positive int AND ``n_bootstrap > 0``. The + (``by_path`` is a positive int OR ``paths_of_interest`` is + set) AND ``n_bootstrap > 0``. The band itself is applied per-horizon as ``cband_conf_int`` on ``path_effects[path]["horizons"][l]`` and rendered as ``cband_lower`` / ``cband_upper`` columns on ``to_dataframe(level="by_path")``. Empty-state contract: - ``None`` when not requested (no bootstrap or ``by_path is None``); - ``{}`` when requested but no path passed both gates (``>=2`` + ``None`` when not requested (no bootstrap, or both ``by_path`` + and ``paths_of_interest`` are ``None``); ``{}`` when requested + but no path passed both gates (``>=2`` valid horizons with finite bootstrap SE ``> 0`` AND a strict majority — more than 50% — of finite sup-t draws). Bands cover joint inference WITHIN a @@ -585,9 +593,9 @@ class ChaisemartinDHaultfoeuilleResults: # conservative upper bound (sum of per-horizon component SEs, # NaN-consistent), matching the global `linear_trends_effects` # convention. - path_cumulated_event_study: Optional[ - Dict[Tuple[int, ...], Dict[int, Dict[str, Any]]] - ] = field(default=None, repr=False) + path_cumulated_event_study: Optional[Dict[Tuple[int, ...], Dict[int, Dict[str, Any]]]] = field( + default=None, repr=False + ) # Per-path joint sup-t simultaneous-band metadata. Keyed by path # tuple; each entry holds `{"crit_value", "alpha", "n_bootstrap", # "method", "n_valid_horizons"}`. Populated when `by_path` is a @@ -1259,13 +1267,26 @@ def _render_path_effects_section( if self.path_effects is None: return if not self.path_effects: + # Distinguish the two empty causes for paths_of_interest + # users (every requested path unobserved) from by_path=k + # users (no panel path has a complete window). + poi = getattr(self._estimator_ref, "paths_of_interest", None) + if poi is not None: + detail_lines = [ + " Every path in paths_of_interest was unobserved or had a window outside L_max+1.", + " (See per-path 'zero observed groups' UserWarnings emitted at fit().)", + ] + else: + detail_lines = [ + " No observed paths have a complete [F_g-1, F_g-1+L_max] window.", + " (See UserWarning emitted at fit(); by_path was a no-op on this panel.)", + ] lines.extend( [ thin, - "Treatment-Path Disaggregation (by_path)".center(width), + "Treatment-Path Disaggregation".center(width), thin, - " No observed paths have a complete [F_g-1, F_g-1+L_max] window.", - " (See UserWarning emitted at fit(); by_path was a no-op " "on this panel.)", + *detail_lines, thin, "", ] @@ -1274,14 +1295,15 @@ def _render_path_effects_section( lines.extend( [ thin, - "Treatment-Path Disaggregation (by_path)".center(width), + "Treatment-Path Disaggregation".center(width), thin, ] ) - for path in sorted( - self.path_effects.keys(), - key=lambda p: self.path_effects[p]["frequency_rank"], - ): + # Iterate in path_effects insertion order so summary preserves + # the user-specified path order under `paths_of_interest`. Under + # `by_path=k`, insertion order matches descending frequency_rank + # (the enumeration sorts by count), so the rendering is identical. + for path in self.path_effects.keys(): entry = self.path_effects[path] rank = entry["frequency_rank"] n_groups = entry["n_groups"] @@ -1337,9 +1359,7 @@ def _render_path_effects_section( ): cum_horizons = self.path_cumulated_event_study[path] if cum_horizons: - lines.append( - " Cumulated Level Effects (DID^{fd}, trends_linear):" - ) + lines.append(" Cumulated Level Effects (DID^{fd}, trends_linear):") for l_h in sorted(cum_horizons.keys()): ce = cum_horizons[l_h] lines.append( @@ -1434,8 +1454,9 @@ def to_dataframe(self, level: str = "overall") -> pd.DataFrame: Available when ``trends_linear=True``. - ``"design2"``: Design-2 switch-in/switch-out descriptive summary. Available when ``design2=True``. - - ``"by_path"``: one row per (path, horizon) when - ``by_path=k`` was passed to the estimator. Columns: + - ``"by_path"``: one row per (path, horizon) when either + ``by_path=k`` or ``paths_of_interest=[(...), ...]`` was + passed to the estimator. Columns: ``path``, ``frequency_rank``, ``n_groups``, ``horizon``, ``effect``, ``se``, ``t_stat``, ``p_value``, ``conf_int_lower``, ``conf_int_upper``, ``n_obs``, @@ -1693,9 +1714,11 @@ def to_dataframe(self, level: str = "overall") -> pd.DataFrame: # Mirrors the linear_trends pattern above. if self.path_effects is None: raise ValueError( - "Path effects not available. Pass by_path=k (positive int) " + "Path effects not available. Pass by_path=k " + "(positive int) or paths_of_interest=[(...), ...] " "to ChaisemartinDHaultfoeuille(drop_larger_lower=False, " - "by_path=k) and L_max >= 1 to fit()." + "by_path=k) (or paths_of_interest=...) and L_max >= 1 " + "to fit()." ) if not self.path_effects: return pd.DataFrame( @@ -1718,10 +1741,11 @@ def to_dataframe(self, level: str = "overall") -> pd.DataFrame: ] ) rows = [] - for path in sorted( - self.path_effects.keys(), - key=lambda p: self.path_effects[p]["frequency_rank"], - ): + # Iterate in path_effects insertion order so the long-format + # table preserves the user-specified path order under + # `paths_of_interest`. Under `by_path=k`, insertion order + # matches descending frequency_rank, so output is identical. + for path in self.path_effects.keys(): entry = self.path_effects[path] rank = entry["frequency_rank"] n_groups = entry["n_groups"] diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 7e66f9b3..f81f7d57 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -242,7 +242,8 @@ ChaisemartinDHaultfoeuille( placebo: bool = True, # Auto-compute single-lag placebo twfe_diagnostic: bool = True, # Auto-compute Theorem 1 TWFE decomposition drop_larger_lower: bool = True, # Drop multi-switch groups (matches R DIDmultiplegtDYN) - by_path: int | None = None, # Top-k per-path event study; requires drop_larger_lower=False, L_max>=1 + by_path: int | None = None, # Top-k per-path event study; requires drop_larger_lower=False, L_max>=1; supports binary or integer-coded discrete D (D in Z); mutex with paths_of_interest + paths_of_interest: list[tuple[int, ...]] | None = None, # User-specified path subset, alternative to by_path=k (Python-only API; mutex with by_path) rank_deficient_action: str = "warn", # Used by TWFE diagnostic OLS ) ``` diff --git a/docs/api/chaisemartin_dhaultfoeuille.rst b/docs/api/chaisemartin_dhaultfoeuille.rst index 576c5c6d..04084618 100644 --- a/docs/api/chaisemartin_dhaultfoeuille.rst +++ b/docs/api/chaisemartin_dhaultfoeuille.rst @@ -19,7 +19,10 @@ integration on placebos; survey support via Taylor-series linearization ``by_path=k`` (mirrors R ``did_multiplegt_dyn(..., by_path=k)``, including per-path backward placebos and per-path joint sup-t simultaneous bands when ``n_bootstrap > 0`` — Python-only extension -beyond R, which provides no joint bands at any surface). +beyond R, which provides no joint bands at any surface) or via +``paths_of_interest=[(...), ...]`` for an explicit user-specified +path subset (Python-only API; mutex with ``by_path``). ``by_path`` +supports binary or integer-coded discrete (D in Z) treatment. The estimator: diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index a59949e6..2b3a29d6 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -638,9 +638,9 @@ The guard is fired by `_survey_se_from_group_if` (analytical and replicate) and - **Note (Phase 3 Design-2 switch-in/switch-out):** Convenience wrapper for Web Appendix Section 1.6 (Assumption 16). Identifies groups with exactly 2 treatment changes (join then leave), reports switch-in and switch-out mean effects. This is a descriptive summary, not a full re-estimation with specialized control pools as described in the paper. **Always uses raw (unadjusted) outcomes** regardless of active `controls`, `trends_linear`, or `trends_nonparam` options - those adjustments apply to the main estimator surface but not to the Design-2 descriptive block. For full adjusted Design-2 estimation with proper control pools, the paper recommends "running the command on a restricted subsample and using `trends_nonparam` for the entry-timing grouping." Activated via `design2=True` in `fit()`, requires `drop_larger_lower=False` to retain 2-switch groups. -- **Note (Phase 3 `by_path` per-path event-study disaggregation):** Per-path disaggregation of the multi-horizon event study, mirroring R `did_multiplegt_dyn(..., by_path=k)`. Activated via `ChaisemartinDHaultfoeuille(by_path=k, drop_larger_lower=False)` where `k` is a positive integer (top-k most common observed paths by switcher-group frequency). **Window convention:** the path tuple for a switcher group `g` is `(D_{g, F_g-1}, D_{g, F_g}, ..., D_{g, F_g-1+L_max})` — length `L_max + 1`, matching R's window `[F_{g-1}, F_{g-1+l}]`. **Ranking:** paths are ranked by descending frequency; ties are broken lexicographically on the path tuple for deterministic ordering, so every selected path has a unique `frequency_rank`. If `by_path` exceeds the number of observed paths, all observed paths are returned with a `UserWarning`. **Per-path SE convention (joiners/leavers precedent):** the per-path influence function follows the joiners-only / leavers-only IF construction at `chaisemartin_dhaultfoeuille.py:5495-5504`: the switcher-side contribution `+S_g * (Y_{g,out} - Y_{g,ref})` is zeroed for groups whose observed trajectory is NOT the selected path; control contributions and the full cohort structure `(D_{g,1}, F_g, S_g)` are unchanged. After applying the singleton-baseline eligible mask and cohort-recentering with the original cohort IDs, the plug-in SE uses the path-specific divisor `N_l_path` (count of path switchers eligible at horizon `l`) — same pattern as `joiners_se` using `joiner_total`. This gives the **within-path mean** estimand `DID_{path,l}` as the within-path average of `DID_{g,l}`. **Degenerate-cohort behavior per path:** when a path's centered IF at some horizon is identically zero (every variance-eligible path switcher forms its own `(D_{g,1}, F_g, S_g)` cohort, or the path has a single contributing group), SE / t_stat / p_value / conf_int are NaN-consistent and a `UserWarning` is emitted scoped to `(path, horizon)`. This mirrors the overall-path degenerate-cohort surface and is common for rare paths with few contributing groups. **Empty-state contract:** `results.path_effects` distinguishes "not requested" (`None`) from "requested but empty" (`{}` — all switchers have windows outside the panel or unobserved cells). The empty-dict case emits a `UserWarning` at fit-time and renders as an explicit "no observed paths" notice in `summary()`; `to_dataframe(level="by_path")` returns an empty DataFrame with the canonical column set (mirrors the `linear_trends` pattern when `trends_linear=True` but no horizons survive). **Requirements:** `drop_larger_lower=False` (multi-switch groups are the object of interest; default `True` filters them out) and `L_max >= 1` (path window depends on the horizon). **Scope:** binary treatment only; combinations with `heterogeneity`, `design2`, `honest_did`, and `survey_design` remain gated behind explicit `NotImplementedError` (deferred to follow-up wave PRs). `n_bootstrap > 0` is now supported — see the **Bootstrap SE** paragraph below. `placebo=True` is now supported per-path — see the **Per-path placebos** paragraph below. **TWFE diagnostic** remains a sample-level summary (not computed per path) in this release. Results are exposed on `results.path_effects` as `Dict[Tuple[int, ...], Dict[str, Any]]` with nested `horizons` dicts per horizon `l`, and on `results.to_dataframe(level="by_path")` as a long-format table with columns `[path, frequency_rank, n_groups, horizon, effect, se, t_stat, p_value, conf_int_lower, conf_int_upper, n_obs, cband_lower, cband_upper, cumulated_effect, cumulated_se]` (the `cband_*` columns are added by the joint sup-t Note below, populated for positive-horizon rows of paths with a finite sup-t crit and NaN otherwise; the `cumulated_*` columns are added by the per-path linear-trends Note below, populated for positive-horizon rows when `trends_linear=True` is set and NaN otherwise). Gated tests live in `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathGates` / `::TestByPathBehavior` / `::TestByPathEdgeCases`. **R-parity** against `DIDmultiplegtDYN 2.3.3` is confirmed at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPath` via two scenarios: `mixed_single_switch_by_path` (2 paths, `by_path=2`) and `multi_path_reversible_by_path` (4 paths, `by_path=3`; path-assignment deterministic on `F_g` so each `(D_{g,1}, F_g, S_g)` cohort contains switchers from a single path). Per-path point estimates and per-path switcher counts match R exactly; per-path SE matches within the Phase 2 multi-horizon SE envelope (observed rtol ≤ 10.2% on the 2-path mixed scenario, ≤ 4.2% on the 4-path cohort-clean scenario). **Deviation from R (cross-path cohort-sharing SE):** our analytical SE is the marginal variance of the path-contribution estimator cohort-centered on the *full-panel* cohort structure (joiners/leavers precedent — non-path switchers contribute to cohort means via their zeroed switcher row). R's `did_multiplegt_dyn(..., by_path=k)` re-runs the estimator per path, so cohort means are computed over the path's own switchers only. When a cohort `(D_{g,1}, F_g, S_g)` spans multiple observed paths, Python and R SE diverge materially (our empirical probes with random post-window toggling saw rtol > 100%); when every cohort is single-path (scenario 13 by design, scenario 14 by construction), the two approaches coincide up to the documented Phase 2 envelope. Practitioners with cohort structures that mix paths should interpret the per-path SE as a within-full-panel marginal variance, not a per-path conditional variance. **Bootstrap SE:** when `n_bootstrap > 0` is set, the top-k paths are enumerated once on the observed data (R-faithful: matches `did_multiplegt_dyn(..., by_path=k, bootstrap=B)`'s path-stability convention — verified empirically against DIDmultiplegtDYN 2.3.3) and the multiplier bootstrap (`bootstrap_weights ∈ {"rademacher", "mammen", "webb"}`) runs per `(path, horizon)` target via the shared `_bootstrap_one_target` / `compute_effect_bootstrap_stats` helpers. Point estimates are unchanged from the analytical path. Bootstrap SE replaces the analytical SE in `path_effects[path]["horizons"][l]["se"]`, and `p_value` / `conf_int` are taken as the **bootstrap percentile** statistics, matching the Round-10 library convention for overall / joiners / leavers / multi-horizon bootstrap (see the `Note (bootstrap inference surface)` elsewhere in this file and the pinned regression `test_bootstrap_p_value_and_ci_propagated_to_top_level`). `t_stat` is SE-derived via `safe_inference` per the anti-pattern rule. Interpretation: inference is *conditional on the observed path set*. **SE inherits the analytical cross-path cohort-sharing deviation:** the bootstrap input is the exact same full-panel cohort-centered path IF that the analytical path computes (`_collect_path_bootstrap_inputs` reuses the same enumeration / cohort IDs / IF construction), so the bootstrap SE is a Monte Carlo analog of the analytical SE — it inherits the same cross-path cohort-sharing deviation from R's per-path re-run convention documented above. On single-path-cohort panels (scenarios 13 and 14 of the R-parity fixture, and any DGP where `(D_{g,1}, F_g, S_g)` cohorts never span multiple observed paths), bootstrap SE tracks analytical SE up to Monte Carlo noise and both coincide with R up to the Phase 2 envelope. On cross-path cohort panels, bootstrap SE inherits the >100% rtol divergence from R that analytical already has. **Deviation from R (CI method):** R's per-path CI is normal-theory around the bootstrap SE (half-width ≈ `1.96·se`); ours is the bootstrap percentile CI, intentionally diverging from R to keep the dCDH inference surface internally consistent across all bootstrap targets. Practitioners who want *unconditional* inference capturing path-selection uncertainty need a pairs-bootstrap (deferred — no R precedent). Positive regressions live in `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathBootstrap` (gated `@pytest.mark.slow`): point-estimate invariance, finite positive SE on non-degenerate panels, SE-within-30%-rtol of analytical on cohort-clean fixtures, degenerate-cohort NaN propagation, Rademacher/Mammen/Webb parity, seed reproducibility, and percentile-vs-normal-theory CI pinning. **Per-path placebos:** when `placebo=True` (and `L_max >= 1`) is combined with `by_path=k`, per-path backward-horizon placebos `DID^{pl}_{path, l}` for `l = 1..L_max` are computed using the same joiners/leavers IF precedent applied to `_compute_per_group_if_placebo_horizon` (with the new `switcher_subset_mask` parameter): switcher contributions are zeroed for groups not in the path; the control pool and the variance-eligible cohort structure `(D_{g,1}, F_g, S_g)` are unchanged. Plug-in SE uses the path-specific divisor `N^{pl}_{l, path}` (count of path switchers eligible at backward lag `l`). Surfaced on `results.path_placebo_event_study[path][-l]` with the same `{effect, se, t_stat, p_value, conf_int, n_obs}` shape as `placebo_event_study` (negative-int inner keys parallel the existing per-path event-study positive-int keys, so a unified forward+backward view is well-formed). **Inherits the cross-path cohort-sharing SE deviation from R** documented above for `path_effects` (same convention applied backward); tracks R within numerical tolerance on single-path-cohort panels and diverges on cohort-mixed panels. Multiplier bootstrap (when `n_bootstrap > 0`) runs per `(path, lag)` target via the same `_bootstrap_one_target` dispatch used for the per-path event-study, with the canonical NaN-on-invalid contract. The bootstrap SE is a Monte Carlo analog of the analytical placebo SE — same per-path centered IF input — and inherits the same deviation. Surfaced through `summary()` (negative-keyed rows rendered alongside positive-keyed event-study rows under each path block) and `to_dataframe(level="by_path")` (`horizon` column takes negative ints for placebo rows). **Empty-state contract:** `results.path_placebo_event_study` mirrors `path_effects` — `None` when `by_path + placebo` was not requested, `{}` when requested but no observed path has a complete window within the panel (same regime that returns `{}` for `path_effects`, with the same fit-time `UserWarning`). R-parity is confirmed at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathPlacebo` on the `multi_path_reversible_by_path_placebo` scenario; positive analytical + bootstrap invariants live in `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathPlacebo` (with the gated `::TestByPathPlacebo::TestBootstrap` subclass). **Per-path covariate residualization (DID^X):** when `controls=[...]` is set with `by_path=k`, the per-baseline OLS residualization (Web Appendix Section 1.2) runs once on the first-differenced outcome BEFORE path enumeration. All four downstream surfaces — analytical per-path SE, bootstrap SE, per-path placebos, and per-path joint sup-t bands — consume the residualized `Y_mat` automatically (Frisch-Waugh-Lovell). Per-period effects remain unadjusted, consistent with the existing `controls` + per-period DID contract (per-period DID does not support residualization). Failed-stratum baselines (rank-deficient X) zero out `N_mat` for affected groups, which the path enumeration treats as ineligible per its existing convention. **Deviation from R on multi-baseline switcher panels (point estimates):** R `did_multiplegt_dyn(..., by_path, controls)` re-runs the per-baseline residualization on each path's restricted subsample (`R/R/did_multiplegt_dyn.R` lines 401-405: rows of the path's switchers OR rows where `yet_to_switch=1 AND baseline matches the path's baseline`). The first-stage residualization sample R uses for path B equals: pre-switch rows of all switchers with matching baseline + all rows of never-switchers with matching baseline — bit-identical to our global first-stage sample under single-baseline switcher panels (every switcher shares the same `D_{g,1}`, regardless of how `F_g` or path identity varies across switchers). Per-path point estimates therefore coincide with R on those panels up to the existing **DID^X first-stage cell-weighting deviation** documented above in `Note (Phase 3 DID^X covariate adjustment)` (Python's first-stage OLS uses equal cell weights — one observation per `(g, t)` cell, consistent with the library's cell-aggregated input convention; R weights by `N_gt`). On panels with one observation per `(g, t)` cell (the common case after the cell-aggregation step in `fit()`), Python matches R bit-exactly: the `multi_path_reversible_by_path_controls` parity fixture has 4 paths with switcher `F_g` values spanning [0..6] under `D_{g,1}=0` and Python matches R to rtol ~1e-11. On multi-baseline switcher panels (some switchers have `D_{g,1}=0`, others have `D_{g,1}=1`) R's per-path subset drops switchers whose baseline differs from the path's baseline, so the per-baseline regression coefficients diverge per path under R and point estimates can diverge between Python and R — a `UserWarning` is emitted at fit-time when this configuration is detected so practitioners do not silently consume estimates that disagree with R. The warning filters to switcher groups only; never-switchers (never-treated + always-treated controls) at multiple baseline values do NOT trigger the warning because they don't affect R's per-path subset construction. **Inherits the cross-path cohort-sharing SE deviation from R** documented above for `path_effects` — bootstrap SE, placebo SE, and sup-t crit are Monte Carlo / joint-distribution analogs of the same residualized analytical IF and carry the same deviation. R-parity is confirmed against `did_multiplegt_dyn(..., by_path=3, controls="X1")` at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathControls` on the `multi_path_reversible_by_path_controls` scenario (single-baseline DGP, exact point-estimate match measured rtol ~1e-11); cross-surface inheritance and the multi-baseline warning are regression-tested at `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathControls` (analytical + bootstrap + placebo + sup-t + `to_dataframe(level="by_path")` cband columns + multi-baseline `UserWarning`). **Per-path linear-trends DID^{fd}:** when `trends_linear=True` is set with `by_path=k`, the first-differencing transform at `chaisemartin_dhaultfoeuille.py:1599-1630` runs once globally BEFORE path enumeration (replaces `Y_mat` with `Z_mat = Y_t - Y_{t-1}` and shrinks the time axis by one), so per-path raw second-differences `DID^{fd}_{path, l}` surface on `path_effects[path]["horizons"][l]` automatically. Per-path cumulated level effects `delta_{path, l} = sum_{l'=1..l} DID^{fd}_{path, l'}` (the quantity R returns under `did_multiplegt_dyn(..., by_path, trends_lin)` per the existing parity test pivot at `tests/test_chaisemartin_dhaultfoeuille_parity.py:403-409`) surface on the new `results.path_cumulated_event_study[path][l]` field — a per-group running sum of `DID^{fd}_{g, l'}` averaged over the path's switchers eligible at horizon `l`, mirroring the global `linear_trends_effects` cumulation logic at `chaisemartin_dhaultfoeuille.py:3340-3398`. SE on the cumulated layer is the conservative upper bound (sum of per-horizon component SEs from `path_effects[path]["horizons"][l]["se"]`, NaN-consistent: any non-finite component yields a NaN cumulated SE). **Post-bootstrap recomputation:** the cumulated layer is built AFTER the bootstrap propagation block at `chaisemartin_dhaultfoeuille.py:3034-3081` so it reads the FINAL post-bootstrap per-horizon SEs (mirrors the global `linear_trends_effects` placement). When `n_bootstrap > 0`, cumulated SE / t / p / CI are derived from bootstrap per-horizon SEs; when bootstrap produces non-finite SE (e.g., `n_bootstrap=1` degenerate distribution), the cumulated layer's full inference tuple is NaN per the library-wide NaN-on-invalid bootstrap contract. `to_dataframe(level="by_path")` exposes `cumulated_effect` and `cumulated_se` columns (always present, NaN-when-None — mirrors the `cband_*` always-present convention from PR #374). `summary()` renders a `Cumulated Level Effects (DID^{fd}, trends_linear)` sub-section under each per-path block. **Path enumeration uses the post-first-differenced `N_mat_fd`**: switchers with `F_g==2` fail the window-eligibility check and are dropped from path enumeration entirely (the existing global `F_g >= 3` warning at line 1620 surfaces the issue), so a path whose switchers all have `F_g < 3` is silently absent from `path_effects` rather than present-with-NaN. **F_g=3 boundary-case divergence (`by_path + trends_linear`):** `F_g=3` switchers have exactly 2 pre-switch periods, which after first-differencing and the `time==1` filter leaves only 1 valid pre-window Z value. R's per-path full-pipeline call handles this single-pre-period regime differently from Python's global-then-disaggregate architecture, producing 30%+ relative divergence on point estimates for paths whose switchers include `F_g=3` (empirically observed on the parity fixture's earlier `F_g=3` variant). A separate `UserWarning` fires at fit-time when the panel includes any `F_g=3` switcher AND `by_path + trends_linear` is set, mirroring the `F_g < 3` exclusion warning. The shipped parity fixture (`single_baseline_multi_path_by_path_trends_lin`) restricts to `F_g >= 4` exclusively to avoid this regime; per-path R parity is asserted only there. **Placebo under `trends_linear` returns RAW per-horizon values** (no per-path placebo cumulation surface) — verified empirically against the existing `joiners_only_trends_lin` parity fixture: R's per-path Placebo_l matches Python's `path_placebo_event_study[path][-l]` (raw) bit-exactly under non-`by_path` trends_lin. **Deviation from R on multi-baseline switcher panels (point estimates):** R `did_multiplegt_dyn(..., by_path, trends_lin)` re-runs the full pipeline (including first-differencing) on each path's restricted subsample, so it operates on different switcher samples per path when switchers have different baseline values `D_{g,1}`. Python first-differences once globally before path enumeration. On single-baseline switcher panels the two architectures coincide; on multi-baseline switcher panels per-path point estimates can diverge — a `UserWarning` is emitted at fit-time when this configuration is detected so practitioners do not silently consume estimates that disagree with R (mirroring the analogous `by_path + controls` warning). Per-path R parity is confirmed against `did_multiplegt_dyn(..., by_path=3, trends_lin=TRUE, placebo=1)` at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathTrendsLinear` on the `single_baseline_multi_path_by_path_trends_lin` scenario (single-baseline + cohort-single-path + `F_g >= 4` DGP designed to eliminate the multi-baseline divergence, the cross-path cohort-sharing deviation, and the F_g=3 boundary case under R's per-path full-pipeline call). Per-path cumulated point estimates match R bit-exactly (rtol ~1e-9) on event horizons under those conditions; cumulated SE_RTOL is widened to `0.20` (vs `0.12` used for non-cumulated by_path parity) because the conservative upper-bound SE compounds the cross-path cohort-sharing deviation under summation. **Placebo parity is intentionally skipped for `trends_linear`**: R's per-path placebo computation re-runs on the path-restricted subsample with different control eligibility than Python's global-then-disaggregate architecture surfaces, producing a sign-and-magnitude divergence on paths whose switchers have minimal pre-window depth (e.g., `F_g=4` switchers). Placebo under `by_path + trends_linear` is exercised via internal regression in `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathTrendsLinear` (finite values, bootstrap inheritance) but not pinned to R bit-by-bit. Cross-surface invariants (analytical + bootstrap + placebo + sup-t + `path_cumulated_event_study` + `to_dataframe` columns + `summary()` rendering) are regression-tested at `TestByPathTrendsLinear`. **Per-path state-set trends:** when `trends_nonparam="state_col"` is set with `by_path=k`, the set membership column is validated and stored once globally as `set_ids_arr` (time-invariance, NaN rejection, partition-coarseness checks unchanged from the non-by_path path). The `set_ids` parameter is threaded through the four per-path IF helpers (`_compute_path_effects`, `_compute_path_placebos`, `_collect_path_bootstrap_inputs`, `_collect_path_placebo_bootstrap_inputs`) so per-path analytical SE, bootstrap, placebos, and sup-t bands all consume the set-restricted control pool automatically. R does NOT first-difference and does NOT cumulate under `trends_nonparam` (unlike `trends_lin`); per-horizon `Effect_l` is a normal DID with set-restricted controls. Per-path R parity is confirmed against `did_multiplegt_dyn(..., by_path=3, trends_nonparam="state", placebo=1)` at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathTrendsNonparam` on the `multi_path_reversible_by_path_trends_nonparam` scenario; per-path point estimates AND placebos match R bit-exactly (rtol ~1e-9), per-path SE matches within the Phase 2 envelope (~13% rtol observed). Cross-surface invariants are regression-tested at `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathTrendsNonparam`. +- **Note (Phase 3 `by_path` per-path event-study disaggregation):** Per-path disaggregation of the multi-horizon event study, mirroring R `did_multiplegt_dyn(..., by_path=k)`. Activated via `ChaisemartinDHaultfoeuille(by_path=k, drop_larger_lower=False)` where `k` is a positive integer (top-k most common observed paths by switcher-group frequency). **Window convention:** the path tuple for a switcher group `g` is `(D_{g, F_g-1}, D_{g, F_g}, ..., D_{g, F_g-1+L_max})` — length `L_max + 1`, matching R's window `[F_{g-1}, F_{g-1+l}]`. **Ranking:** paths are ranked by descending frequency; ties are broken lexicographically on the path tuple for deterministic ordering, so every selected path has a unique `frequency_rank`. If `by_path` exceeds the number of observed paths, all observed paths are returned with a `UserWarning`. **Per-path SE convention (joiners/leavers precedent):** the per-path influence function follows the joiners-only / leavers-only IF construction at `chaisemartin_dhaultfoeuille.py:5495-5504`: the switcher-side contribution `+S_g * (Y_{g,out} - Y_{g,ref})` is zeroed for groups whose observed trajectory is NOT the selected path; control contributions and the full cohort structure `(D_{g,1}, F_g, S_g)` are unchanged. After applying the singleton-baseline eligible mask and cohort-recentering with the original cohort IDs, the plug-in SE uses the path-specific divisor `N_l_path` (count of path switchers eligible at horizon `l`) — same pattern as `joiners_se` using `joiner_total`. This gives the **within-path mean** estimand `DID_{path,l}` as the within-path average of `DID_{g,l}`. **Degenerate-cohort behavior per path:** when a path's centered IF at some horizon is identically zero (every variance-eligible path switcher forms its own `(D_{g,1}, F_g, S_g)` cohort, or the path has a single contributing group), SE / t_stat / p_value / conf_int are NaN-consistent and a `UserWarning` is emitted scoped to `(path, horizon)`. This mirrors the overall-path degenerate-cohort surface and is common for rare paths with few contributing groups. **Empty-state contract:** `results.path_effects` distinguishes "not requested" (`None`) from "requested but empty" (`{}` — all switchers have windows outside the panel or unobserved cells). The empty-dict case emits a `UserWarning` at fit-time and renders as an explicit "no observed paths" notice in `summary()`; `to_dataframe(level="by_path")` returns an empty DataFrame with the canonical column set (mirrors the `linear_trends` pattern when `trends_linear=True` but no horizons survive). **Requirements:** `drop_larger_lower=False` (multi-switch groups are the object of interest; default `True` filters them out) and `L_max >= 1` (path window depends on the horizon). **Scope:** combinations with `heterogeneity`, `design2`, `honest_did`, and `survey_design` remain gated behind explicit `NotImplementedError` (deferred to follow-up wave PRs). `n_bootstrap > 0` is now supported — see the **Bootstrap SE** paragraph below. `placebo=True` is now supported per-path — see the **Per-path placebos** paragraph below. **TWFE diagnostic** remains a sample-level summary (not computed per path) in this release. Results are exposed on `results.path_effects` as `Dict[Tuple[int, ...], Dict[str, Any]]` with nested `horizons` dicts per horizon `l`, and on `results.to_dataframe(level="by_path")` as a long-format table with columns `[path, frequency_rank, n_groups, horizon, effect, se, t_stat, p_value, conf_int_lower, conf_int_upper, n_obs, cband_lower, cband_upper, cumulated_effect, cumulated_se]` (the `cband_*` columns are added by the joint sup-t Note below, populated for positive-horizon rows of paths with a finite sup-t crit and NaN otherwise; the `cumulated_*` columns are added by the per-path linear-trends Note below, populated for positive-horizon rows when `trends_linear=True` is set and NaN otherwise). Gated tests live in `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathGates` / `::TestByPathBehavior` / `::TestByPathEdgeCases`. **R-parity** against `DIDmultiplegtDYN 2.3.3` is confirmed at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPath` via two scenarios: `mixed_single_switch_by_path` (2 paths, `by_path=2`) and `multi_path_reversible_by_path` (4 paths, `by_path=3`; path-assignment deterministic on `F_g` so each `(D_{g,1}, F_g, S_g)` cohort contains switchers from a single path). Per-path point estimates and per-path switcher counts match R exactly; per-path SE matches within the Phase 2 multi-horizon SE envelope (observed rtol ≤ 10.2% on the 2-path mixed scenario, ≤ 4.2% on the 4-path cohort-clean scenario). **Deviation from R (cross-path cohort-sharing SE):** our analytical SE is the marginal variance of the path-contribution estimator cohort-centered on the *full-panel* cohort structure (joiners/leavers precedent — non-path switchers contribute to cohort means via their zeroed switcher row). R's `did_multiplegt_dyn(..., by_path=k)` re-runs the estimator per path, so cohort means are computed over the path's own switchers only. When a cohort `(D_{g,1}, F_g, S_g)` spans multiple observed paths, Python and R SE diverge materially (our empirical probes with random post-window toggling saw rtol > 100%); when every cohort is single-path (scenario 13 by design, scenario 14 by construction), the two approaches coincide up to the documented Phase 2 envelope. Practitioners with cohort structures that mix paths should interpret the per-path SE as a within-full-panel marginal variance, not a per-path conditional variance. **Bootstrap SE:** when `n_bootstrap > 0` is set, the top-k paths are enumerated once on the observed data (R-faithful: matches `did_multiplegt_dyn(..., by_path=k, bootstrap=B)`'s path-stability convention — verified empirically against DIDmultiplegtDYN 2.3.3) and the multiplier bootstrap (`bootstrap_weights ∈ {"rademacher", "mammen", "webb"}`) runs per `(path, horizon)` target via the shared `_bootstrap_one_target` / `compute_effect_bootstrap_stats` helpers. Point estimates are unchanged from the analytical path. Bootstrap SE replaces the analytical SE in `path_effects[path]["horizons"][l]["se"]`, and `p_value` / `conf_int` are taken as the **bootstrap percentile** statistics, matching the Round-10 library convention for overall / joiners / leavers / multi-horizon bootstrap (see the `Note (bootstrap inference surface)` elsewhere in this file and the pinned regression `test_bootstrap_p_value_and_ci_propagated_to_top_level`). `t_stat` is SE-derived via `safe_inference` per the anti-pattern rule. Interpretation: inference is *conditional on the observed path set*. **SE inherits the analytical cross-path cohort-sharing deviation:** the bootstrap input is the exact same full-panel cohort-centered path IF that the analytical path computes (`_collect_path_bootstrap_inputs` reuses the same enumeration / cohort IDs / IF construction), so the bootstrap SE is a Monte Carlo analog of the analytical SE — it inherits the same cross-path cohort-sharing deviation from R's per-path re-run convention documented above. On single-path-cohort panels (scenarios 13 and 14 of the R-parity fixture, and any DGP where `(D_{g,1}, F_g, S_g)` cohorts never span multiple observed paths), bootstrap SE tracks analytical SE up to Monte Carlo noise and both coincide with R up to the Phase 2 envelope. On cross-path cohort panels, bootstrap SE inherits the >100% rtol divergence from R that analytical already has. **Deviation from R (CI method):** R's per-path CI is normal-theory around the bootstrap SE (half-width ≈ `1.96·se`); ours is the bootstrap percentile CI, intentionally diverging from R to keep the dCDH inference surface internally consistent across all bootstrap targets. Practitioners who want *unconditional* inference capturing path-selection uncertainty need a pairs-bootstrap (deferred — no R precedent). Positive regressions live in `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathBootstrap` (gated `@pytest.mark.slow`): point-estimate invariance, finite positive SE on non-degenerate panels, SE-within-30%-rtol of analytical on cohort-clean fixtures, degenerate-cohort NaN propagation, Rademacher/Mammen/Webb parity, seed reproducibility, and percentile-vs-normal-theory CI pinning. **Per-path placebos:** when `placebo=True` (and `L_max >= 1`) is combined with `by_path=k`, per-path backward-horizon placebos `DID^{pl}_{path, l}` for `l = 1..L_max` are computed using the same joiners/leavers IF precedent applied to `_compute_per_group_if_placebo_horizon` (with the new `switcher_subset_mask` parameter): switcher contributions are zeroed for groups not in the path; the control pool and the variance-eligible cohort structure `(D_{g,1}, F_g, S_g)` are unchanged. Plug-in SE uses the path-specific divisor `N^{pl}_{l, path}` (count of path switchers eligible at backward lag `l`). Surfaced on `results.path_placebo_event_study[path][-l]` with the same `{effect, se, t_stat, p_value, conf_int, n_obs}` shape as `placebo_event_study` (negative-int inner keys parallel the existing per-path event-study positive-int keys, so a unified forward+backward view is well-formed). **Inherits the cross-path cohort-sharing SE deviation from R** documented above for `path_effects` (same convention applied backward); tracks R within numerical tolerance on single-path-cohort panels and diverges on cohort-mixed panels. Multiplier bootstrap (when `n_bootstrap > 0`) runs per `(path, lag)` target via the same `_bootstrap_one_target` dispatch used for the per-path event-study, with the canonical NaN-on-invalid contract. The bootstrap SE is a Monte Carlo analog of the analytical placebo SE — same per-path centered IF input — and inherits the same deviation. Surfaced through `summary()` (negative-keyed rows rendered alongside positive-keyed event-study rows under each path block) and `to_dataframe(level="by_path")` (`horizon` column takes negative ints for placebo rows). **Empty-state contract:** `results.path_placebo_event_study` mirrors `path_effects` — `None` when `by_path + placebo` was not requested, `{}` when requested but no observed path has a complete window within the panel (same regime that returns `{}` for `path_effects`, with the same fit-time `UserWarning`). R-parity is confirmed at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathPlacebo` on the `multi_path_reversible_by_path_placebo` scenario; positive analytical + bootstrap invariants live in `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathPlacebo` (with the gated `::TestByPathPlacebo::TestBootstrap` subclass). **Per-path covariate residualization (DID^X):** when `controls=[...]` is set with `by_path=k`, the per-baseline OLS residualization (Web Appendix Section 1.2) runs once on the first-differenced outcome BEFORE path enumeration. All four downstream surfaces — analytical per-path SE, bootstrap SE, per-path placebos, and per-path joint sup-t bands — consume the residualized `Y_mat` automatically (Frisch-Waugh-Lovell). Per-period effects remain unadjusted, consistent with the existing `controls` + per-period DID contract (per-period DID does not support residualization). Failed-stratum baselines (rank-deficient X) zero out `N_mat` for affected groups, which the path enumeration treats as ineligible per its existing convention. **Deviation from R on multi-baseline switcher panels (point estimates):** R `did_multiplegt_dyn(..., by_path, controls)` re-runs the per-baseline residualization on each path's restricted subsample (`R/R/did_multiplegt_dyn.R` lines 401-405: rows of the path's switchers OR rows where `yet_to_switch=1 AND baseline matches the path's baseline`). The first-stage residualization sample R uses for path B equals: pre-switch rows of all switchers with matching baseline + all rows of never-switchers with matching baseline — bit-identical to our global first-stage sample under single-baseline switcher panels (every switcher shares the same `D_{g,1}`, regardless of how `F_g` or path identity varies across switchers). Per-path point estimates therefore coincide with R on those panels up to the existing **DID^X first-stage cell-weighting deviation** documented above in `Note (Phase 3 DID^X covariate adjustment)` (Python's first-stage OLS uses equal cell weights — one observation per `(g, t)` cell, consistent with the library's cell-aggregated input convention; R weights by `N_gt`). On panels with one observation per `(g, t)` cell (the common case after the cell-aggregation step in `fit()`), Python matches R bit-exactly: the `multi_path_reversible_by_path_controls` parity fixture has 4 paths with switcher `F_g` values spanning [0..6] under `D_{g,1}=0` and Python matches R to rtol ~1e-11. On multi-baseline switcher panels (some switchers have `D_{g,1}=0`, others have `D_{g,1}=1`) R's per-path subset drops switchers whose baseline differs from the path's baseline, so the per-baseline regression coefficients diverge per path under R and point estimates can diverge between Python and R — a `UserWarning` is emitted at fit-time when this configuration is detected so practitioners do not silently consume estimates that disagree with R. The warning filters to switcher groups only; never-switchers (never-treated + always-treated controls) at multiple baseline values do NOT trigger the warning because they don't affect R's per-path subset construction. **Inherits the cross-path cohort-sharing SE deviation from R** documented above for `path_effects` — bootstrap SE, placebo SE, and sup-t crit are Monte Carlo / joint-distribution analogs of the same residualized analytical IF and carry the same deviation. R-parity is confirmed against `did_multiplegt_dyn(..., by_path=3, controls="X1")` at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathControls` on the `multi_path_reversible_by_path_controls` scenario (single-baseline DGP, exact point-estimate match measured rtol ~1e-11); cross-surface inheritance and the multi-baseline warning are regression-tested at `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathControls` (analytical + bootstrap + placebo + sup-t + `to_dataframe(level="by_path")` cband columns + multi-baseline `UserWarning`). **Per-path linear-trends DID^{fd}:** when `trends_linear=True` is set with `by_path=k`, the first-differencing transform at `chaisemartin_dhaultfoeuille.py:1599-1630` runs once globally BEFORE path enumeration (replaces `Y_mat` with `Z_mat = Y_t - Y_{t-1}` and shrinks the time axis by one), so per-path raw second-differences `DID^{fd}_{path, l}` surface on `path_effects[path]["horizons"][l]` automatically. Per-path cumulated level effects `delta_{path, l} = sum_{l'=1..l} DID^{fd}_{path, l'}` (the quantity R returns under `did_multiplegt_dyn(..., by_path, trends_lin)` per the existing parity test pivot at `tests/test_chaisemartin_dhaultfoeuille_parity.py:403-409`) surface on the new `results.path_cumulated_event_study[path][l]` field — a per-group running sum of `DID^{fd}_{g, l'}` averaged over the path's switchers eligible at horizon `l`, mirroring the global `linear_trends_effects` cumulation logic at `chaisemartin_dhaultfoeuille.py:3340-3398`. SE on the cumulated layer is the conservative upper bound (sum of per-horizon component SEs from `path_effects[path]["horizons"][l]["se"]`, NaN-consistent: any non-finite component yields a NaN cumulated SE). **Post-bootstrap recomputation:** the cumulated layer is built AFTER the bootstrap propagation block at `chaisemartin_dhaultfoeuille.py:3034-3081` so it reads the FINAL post-bootstrap per-horizon SEs (mirrors the global `linear_trends_effects` placement). When `n_bootstrap > 0`, cumulated SE / t / p / CI are derived from bootstrap per-horizon SEs; when bootstrap produces non-finite SE (e.g., `n_bootstrap=1` degenerate distribution), the cumulated layer's full inference tuple is NaN per the library-wide NaN-on-invalid bootstrap contract. `to_dataframe(level="by_path")` exposes `cumulated_effect` and `cumulated_se` columns (always present, NaN-when-None — mirrors the `cband_*` always-present convention from PR #374). `summary()` renders a `Cumulated Level Effects (DID^{fd}, trends_linear)` sub-section under each per-path block. **Path enumeration uses the post-first-differenced `N_mat_fd`**: switchers with `F_g==2` fail the window-eligibility check and are dropped from path enumeration entirely (the existing global `F_g >= 3` warning at line 1620 surfaces the issue), so a path whose switchers all have `F_g < 3` is silently absent from `path_effects` rather than present-with-NaN. **F_g=3 boundary-case divergence (`by_path + trends_linear`):** `F_g=3` switchers have exactly 2 pre-switch periods, which after first-differencing and the `time==1` filter leaves only 1 valid pre-window Z value. R's per-path full-pipeline call handles this single-pre-period regime differently from Python's global-then-disaggregate architecture, producing 30%+ relative divergence on point estimates for paths whose switchers include `F_g=3` (empirically observed on the parity fixture's earlier `F_g=3` variant). A separate `UserWarning` fires at fit-time when the panel includes any `F_g=3` switcher AND `by_path + trends_linear` is set, mirroring the `F_g < 3` exclusion warning. The shipped parity fixture (`single_baseline_multi_path_by_path_trends_lin`) restricts to `F_g >= 4` exclusively to avoid this regime; per-path R parity is asserted only there. **Placebo under `trends_linear` returns RAW per-horizon values** (no per-path placebo cumulation surface) — verified empirically against the existing `joiners_only_trends_lin` parity fixture: R's per-path Placebo_l matches Python's `path_placebo_event_study[path][-l]` (raw) bit-exactly under non-`by_path` trends_lin. **Deviation from R on multi-baseline switcher panels (point estimates):** R `did_multiplegt_dyn(..., by_path, trends_lin)` re-runs the full pipeline (including first-differencing) on each path's restricted subsample, so it operates on different switcher samples per path when switchers have different baseline values `D_{g,1}`. Python first-differences once globally before path enumeration. On single-baseline switcher panels the two architectures coincide; on multi-baseline switcher panels per-path point estimates can diverge — a `UserWarning` is emitted at fit-time when this configuration is detected so practitioners do not silently consume estimates that disagree with R (mirroring the analogous `by_path + controls` warning). Per-path R parity is confirmed against `did_multiplegt_dyn(..., by_path=3, trends_lin=TRUE, placebo=1)` at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathTrendsLinear` on the `single_baseline_multi_path_by_path_trends_lin` scenario (single-baseline + cohort-single-path + `F_g >= 4` DGP designed to eliminate the multi-baseline divergence, the cross-path cohort-sharing deviation, and the F_g=3 boundary case under R's per-path full-pipeline call). Per-path cumulated point estimates match R bit-exactly (rtol ~1e-9) on event horizons under those conditions; cumulated SE_RTOL is widened to `0.20` (vs `0.12` used for non-cumulated by_path parity) because the conservative upper-bound SE compounds the cross-path cohort-sharing deviation under summation. **Placebo parity is intentionally skipped for `trends_linear`**: R's per-path placebo computation re-runs on the path-restricted subsample with different control eligibility than Python's global-then-disaggregate architecture surfaces, producing a sign-and-magnitude divergence on paths whose switchers have minimal pre-window depth (e.g., `F_g=4` switchers). Placebo under `by_path + trends_linear` is exercised via internal regression in `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathTrendsLinear` (finite values, bootstrap inheritance) but not pinned to R bit-by-bit. Cross-surface invariants (analytical + bootstrap + placebo + sup-t + `path_cumulated_event_study` + `to_dataframe` columns + `summary()` rendering) are regression-tested at `TestByPathTrendsLinear`. **Per-path state-set trends:** when `trends_nonparam="state_col"` is set with `by_path=k`, the set membership column is validated and stored once globally as `set_ids_arr` (time-invariance, NaN rejection, partition-coarseness checks unchanged from the non-by_path path). The `set_ids` parameter is threaded through the four per-path IF helpers (`_compute_path_effects`, `_compute_path_placebos`, `_collect_path_bootstrap_inputs`, `_collect_path_placebo_bootstrap_inputs`) so per-path analytical SE, bootstrap, placebos, and sup-t bands all consume the set-restricted control pool automatically. R does NOT first-difference and does NOT cumulate under `trends_nonparam` (unlike `trends_lin`); per-horizon `Effect_l` is a normal DID with set-restricted controls. Per-path R parity is confirmed against `did_multiplegt_dyn(..., by_path=3, trends_nonparam="state", placebo=1)` at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathTrendsNonparam` on the `multi_path_reversible_by_path_trends_nonparam` scenario; per-path point estimates AND placebos match R bit-exactly (rtol ~1e-9), per-path SE matches within the Phase 2 envelope (~13% rtol observed). Cross-surface invariants are regression-tested at `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathTrendsNonparam`. **Per-path non-binary treatment:** integer-coded discrete treatment (D in Z, e.g. ordinal {0, 1, 2}) is supported under `by_path=k` and `paths_of_interest`. Path tuples become integer-state tuples (`(0, 2, 2, 2)`) keyed bit-for-bit against R's comma-separated path strings (`"0,2,2,2"`) for D in {0..9}. Continuous D (e.g. `1.5`) raises `ValueError` at fit-time per the no-silent-failures contract — the existing `int(round(float(v)))` cast in `_enumerate_treatment_paths` is now defensive (no-op for integer-coded D). **Deviation from R for D >= 10:** R's `did_multiplegt_by_path` derives the per-path baseline via `path_index$baseline_XX <- substr(path_index$path, 1, 1)` (extracted 2026-05-03 via `Rscript -e 'cat(paste(deparse(DIDmultiplegtDYN:::did_multiplegt_by_path), collapse="\n"))'`), capturing only the first character of the comma-separated path string. For D >= 10 this captures `"1"` instead of `"12"` for `path = "12,12,..."`, mis-allocating R's per-path control-pool subset. Python's tuple-key matching is correct in this regime; the per-path point estimates we compute are correct, R's per-path subset for the same path is buggy. The shipped parity scenario stays in `D in {0, 1, 2}` to avoid the R bug; R-parity for D in {0..9} is asserted at `tests/test_chaisemartin_dhaultfoeuille_parity.py::TestDCDHDynRParityByPathNonBinary` on the `multi_path_reversible_by_path_non_binary` scenario (78 switchers, 3 paths, single-baseline custom DGP, F_g >= 4) — per-path point estimates match R bit-exactly (rtol ~1e-9 events; rtol+atol envelope for placebo near-zero values), SE inherits the documented cross-path cohort-sharing deviation (~5% rtol observed; SE_RTOL=0.15 envelope). Cross-surface invariants regression-tested at `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathNonBinary`. **Per-path user-specified path selection (`paths_of_interest`):** Python-only API extension — R's `did_multiplegt_dyn(..., by_path=k)` only accepts a positive int (top-k automatic ranking) or `-1` (all observed paths) and provides no list-based selection. Activated via `ChaisemartinDHaultfoeuille(paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)], drop_larger_lower=False)` as an alternative to `by_path=k`; the two are **mutually exclusive** (setting both raises `ValueError` at `__init__` and `set_params` time). Each path tuple must have length `L_max + 1`; the type / element / non-empty / length-uniformity checks fire at `__init__`, the length-vs-L_max check fires at fit-time. `bool` and `np.bool_` are explicitly rejected; `np.integer` is accepted and canonicalized to Python `int` for tuple-key consistency. Duplicates emit a `UserWarning` and are deduplicated; paths not observed in the panel emit a `UserWarning` and are omitted from `path_effects`. Paths appear in `results.path_effects` in the user-specified order, modulo deduplication and unobserved-path filtering. Composes with non-binary D and all downstream `by_path` surfaces (bootstrap, per-path placebos, per-path joint sup-t bands, `controls`, `trends_linear`, `trends_nonparam`) — mechanical filter on observed paths, no methodology change. Behavior + cross-feature regressions live at `tests/test_chaisemartin_dhaultfoeuille.py::TestPathsOfInterest`. -- **Note (Phase 3 `by_path` per-path joint sup-t bands):** When `n_bootstrap > 0` is set with `by_path=k`, per-path joint sup-t simultaneous confidence bands are computed across horizons `1..L_max` within each path. **Methodology:** a single `(n_bootstrap, n_eligible)` multiplier weight matrix (using the estimator's configured `bootstrap_weights` — Rademacher / Mammen / Webb) is drawn per path and broadcast across all horizons of that path, producing correlated bootstrap distributions across horizons within the path. The path-specific critical value `c_p = quantile(max_l |t_l|, 1 - α)` is then used to construct symmetric joint bands `effect_l ± c_p · se_l` per horizon, surfaced in `path_effects[path]["horizons"][l]["cband_conf_int"]` and at top-level `results.path_sup_t_bands[path] = {"crit_value", "alpha", "n_bootstrap", "method", "n_valid_horizons"}`. **Gates:** a path must have `>= 2` valid horizons (finite bootstrap SE > 0) AND a strict majority (more than 50%) of finite sup-t draws to receive a band; otherwise the path is absent from `path_sup_t_bands`. Both gates mirror the OVERALL `event_study_sup_t_bands` semantics at `chaisemartin_dhaultfoeuille_bootstrap.py:605,612`: `len(valid_horizons) >= 2` AND `finite_mask.sum() > 0.5 * n_bootstrap`. Exactly half-finite draws are NOT enough — the gate is strictly greater than half. **Empty-state contract:** `path_sup_t_bands is None` when not requested (no bootstrap or `by_path is None`); `{}` when requested but no path passes both gates. **`to_dataframe(level="by_path")` integration:** the table now includes `cband_lower` / `cband_upper` columns for parity with OVERALL `level="event_study"`; populated for positive-horizon rows of paths with a finite sup-t crit, NaN for placebo rows / unbanded paths / the requested-but-empty fallback DataFrame. **Methodology asymmetry vs OVERALL:** OVERALL sup-t reuses the same multi-horizon shared-draw distribution for both the SE in the t-stat denominator and the bootstrap distribution in the numerator. The per-path sup-t draws a fresh shared weight matrix per path AFTER the per-path SE bootstrap block has already populated `results.path_ses` via independent per-(path, horizon) draws — numerator: fresh shared draws, denominator: bootstrap SEs from the earlier independent draws. Asymptotically equivalent to OVERALL's self-consistent reuse, but NOT bit-identical. The fresh draw is intentional: it preserves RNG-state isolation and keeps every existing per-path SE seed-reproducibility test bit-stable post-implementation. **Inherited deviation from R:** the bootstrap SE used as the t-stat denominator carries the cross-path cohort-sharing SE deviation from R documented for `path_effects` above; the per-path sup-t crit therefore inherits the same deviation. **Interpretation:** the band covers joint inference *within a single path across horizons*; it does NOT provide simultaneous coverage *across paths* (a different inference target requiring a `path × horizon` re-derivation, deferred to a future wave). **Deviation from R:** `did_multiplegt_dyn` provides no joint / sup-t / simultaneous bands at any surface — this is a Python-only methodology extension, consistent with the existing OVERALL `event_study_sup_t_bands` (also Python-only). Regression test anchor: `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathSupTBands`. +- **Note (Phase 3 `by_path` per-path joint sup-t bands):** When `n_bootstrap > 0` is set with `by_path=k`, per-path joint sup-t simultaneous confidence bands are computed across horizons `1..L_max` within each path. **Methodology:** a single `(n_bootstrap, n_eligible)` multiplier weight matrix (using the estimator's configured `bootstrap_weights` — Rademacher / Mammen / Webb) is drawn per path and broadcast across all horizons of that path, producing correlated bootstrap distributions across horizons within the path. The path-specific critical value `c_p = quantile(max_l |t_l|, 1 - α)` is then used to construct symmetric joint bands `effect_l ± c_p · se_l` per horizon, surfaced in `path_effects[path]["horizons"][l]["cband_conf_int"]` and at top-level `results.path_sup_t_bands[path] = {"crit_value", "alpha", "n_bootstrap", "method", "n_valid_horizons"}`. **Gates:** a path must have `>= 2` valid horizons (finite bootstrap SE > 0) AND a strict majority (more than 50%) of finite sup-t draws to receive a band; otherwise the path is absent from `path_sup_t_bands`. Both gates mirror the OVERALL `event_study_sup_t_bands` semantics at `chaisemartin_dhaultfoeuille_bootstrap.py:605,612`: `len(valid_horizons) >= 2` AND `finite_mask.sum() > 0.5 * n_bootstrap`. Exactly half-finite draws are NOT enough — the gate is strictly greater than half. **Empty-state contract:** `path_sup_t_bands is None` when not requested (no bootstrap, or both `by_path` and `paths_of_interest` are `None`); `{}` when requested but no path passes both gates. **`to_dataframe(level="by_path")` integration:** the table now includes `cband_lower` / `cband_upper` columns for parity with OVERALL `level="event_study"`; populated for positive-horizon rows of paths with a finite sup-t crit, NaN for placebo rows / unbanded paths / the requested-but-empty fallback DataFrame. **Methodology asymmetry vs OVERALL:** OVERALL sup-t reuses the same multi-horizon shared-draw distribution for both the SE in the t-stat denominator and the bootstrap distribution in the numerator. The per-path sup-t draws a fresh shared weight matrix per path AFTER the per-path SE bootstrap block has already populated `results.path_ses` via independent per-(path, horizon) draws — numerator: fresh shared draws, denominator: bootstrap SEs from the earlier independent draws. Asymptotically equivalent to OVERALL's self-consistent reuse, but NOT bit-identical. The fresh draw is intentional: it preserves RNG-state isolation and keeps every existing per-path SE seed-reproducibility test bit-stable post-implementation. **Inherited deviation from R:** the bootstrap SE used as the t-stat denominator carries the cross-path cohort-sharing SE deviation from R documented for `path_effects` above; the per-path sup-t crit therefore inherits the same deviation. **Interpretation:** the band covers joint inference *within a single path across horizons*; it does NOT provide simultaneous coverage *across paths* (a different inference target requiring a `path × horizon` re-derivation, deferred to a future wave). **Deviation from R:** `did_multiplegt_dyn` provides no joint / sup-t / simultaneous bands at any surface — this is a Python-only methodology extension, consistent with the existing OVERALL `event_study_sup_t_bands` (also Python-only). Regression test anchor: `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathSupTBands`. **Reference implementation(s):** - R: [`DIDmultiplegtDYN`](https://cran.r-project.org/package=DIDmultiplegtDYN) (CRAN, maintained by the paper authors). The Python implementation matches `did_multiplegt_dyn(..., effects=1)` at horizon `l = 1`. Parity tests live in `tests/test_chaisemartin_dhaultfoeuille_parity.py`. @@ -668,7 +668,7 @@ The guard is fired by `_survey_se_from_group_if` (analytical and replicate) and - [x] State-set-specific trends via control-pool restriction (Web Appendix Section 1.4) - [x] Heterogeneity testing via saturated OLS (Web Appendix Section 1.5, Lemma 7) - [x] Design-2 switch-in/switch-out descriptive wrapper (Web Appendix Section 1.6) -- [x] `by_path` per-path event-study disaggregation (binary treatment, joiners/leavers IF precedent; mirrors R `did_multiplegt_dyn(..., by_path=k)`) +- [x] `by_path` per-path event-study disaggregation (binary or integer-coded discrete treatment, joiners/leavers IF precedent; mirrors R `did_multiplegt_dyn(..., by_path=k)`); plus `paths_of_interest=[(...), ...]` for user-specified path subsets (Python-only API; mutex with `by_path`) - [x] HonestDiD (Rambachan-Roth 2023) integration on placebo + event study surface - [x] Survey design support: pweight with strata/PSU/FPC via Taylor Series Linearization (analytical) **or replicate-weight variance (BRR/Fay/JK1/JKn/SDR)**, covering the main ATT surface, covariate adjustment (DID^X), heterogeneity testing, the TWFE diagnostic (fit and standalone `twowayfeweights()` helper), and HonestDiD bounds. Opt-in **PSU-level Hall-Mammen wild bootstrap** is also supported via `n_bootstrap > 0`. - **Note (Survey IF expansion — library convention):** Survey IF expansion is a library extension not in the dCDH papers (the paper's plug-in variance assumes iid sampling). The library convention builds observation-level `psi_i` by proportionally distributing per-group IF mass within weight share: either at the group level (`psi_i = U_centered[g] * w_i / W_g`, the previous convention) or at the per-`(g, t)` cell level via the cell-period allocator shipped in this release. Cell-level expansion: decompose `U[g]` into per-period attributions `U[g, t]`, cohort-center each column independently, then expand to observation level as `psi_i = U_centered_per_period[g_i, t_i] * (w_i / W_{g_i, t_i})`. Binder (1983) stratified-PSU variance aggregates the resulting `psi` at PSU level. **Post-period attribution convention:** each transition term in the IF sum (of the form `role_weight * (Y_{g, t} - Y_{g, t-1})` for DID_M or `S_g * (Y_{g, out} - Y_{g, ref})` for DID_l) is attributed as a single *difference* to the POST-period cell, not split into a `+Y_post` / `-Y_pre` pair across two cells. This is a library *convention*, not a theorem — adopted because it preserves the group-sum, PSU-sum, and cohort-sum identities of the previous group-level expansion (so Binder variance coincides with the group-level variance under the auto-injected `psu=group`) and because Monte Carlo coverage at nominal 95% is empirically close to nominal on a DGP where PSUs vary across the cells of each group (see `tests/test_dcdh_cell_period_coverage.py`). A covariance-aware two-cell allocator is a plausible alternative and may be worth exploring if future designs motivate an explicit observation-level IF derivation; the method currently in the library is **not derived from the observation-level survey linearization of the contrast** and makes no stronger claim than "coverage is approximately nominal under the tested DGPs and the group-sum identity holds exactly." Under within-group-constant PSU (the pre-allocator accepted input), per-cell sums telescope to `U_centered[g]` and Binder variance is byte-identical (up to single-ULP floating-point noise) to the previous group-level expansion. **Strata and PSU must be constant within each `(g, t)` cell** (trivially satisfied in one-obs-per-cell panels — the canonical dCDH structure); variation **across cells of a group** is supported by the allocator. Within-group-varying **weights** are supported as before. When `survey_design.psu` is not specified, `fit()` auto-injects `psu=` so the TSL variance, `df_survey`, and t-based inference match the per-group PSU structure. **Strata that vary across cells of a group require either an explicit `psu=` or the original `SurveyDesign(..., nest=True)` flag** — under `nest=True` the resolver combines `(stratum, psu)` into globally-unique labels, so the auto-injected `psu=` is re-labeled per stratum and the cell allocator proceeds. Only the `nest=False` + varying-strata + omitted-psu combination is rejected up front with a targeted `ValueError` at `fit()` time (the synthesized PSU column would reuse group labels across strata and trip the cross-stratum PSU uniqueness check in `SurveyDesign.resolve()`). Under replicate-weight designs, the same cell-level `psi_i` is aggregated via Rao-Wu weight-ratio rescaling (`compute_replicate_if_variance` at `diff_diff/survey.py:1681`) rather than the Binder TSL formula. All five methods (BRR/Fay/JK1/JKn/SDR) are supported method-agnostically through the unified helper; the effective `df_survey` is reduced to `min(n_valid) - 1` across IF sites when some replicate solves fail (matching `efficient_did.py:1133-1135` and `triple_diff.py:676-686` precedents). Under DID^X, the first-stage residualization coefficient `theta_hat` is computed once on full-sample weights and treated as fixed (FWL plug-in IF convention) — per-replicate refits of `theta_hat` are not performed. **Post-period attribution extends to heterogeneity (Binder TSL branch only):** the heterogeneity WLS coefficient IF `ψ_g = inv(X'WX)[1,:] @ x_g * W_g * r_g` is attributed in full to the single post-period cell `(g, out_idx)` at each horizon (same single-cell convention as DID_l), then expanded as `ψ_i = ψ_g * (w_i / W_{g, out_idx})`, and fed through `compute_survey_if_variance`. Under PSU=group the PSU-level aggregate telescopes to `ψ_g`, so Binder variance is byte-identical relative to the pre-cell-period release; under within-group-varying PSU mass lands in the post-period PSU. **Replicate-weight branch keeps the legacy group-level allocator** `ψ_i = ψ_g * (w_i / W_g)` because `compute_replicate_if_variance` computes `θ_r = sum_i ratio_ir * ψ_i` at observation level and is therefore not PSU-telescoping: redistributing mass onto the post-period cell would silently change the replicate SE whenever a replicate column's ratios vary within a group (the library accepts arbitrary per-row replicate matrices, not just PSU-aligned ones). The legacy allocator preserves byte-identity of the replicate SE for every previously-supported fit. Replicate + within-group-varying PSU is unreachable by construction (`SurveyDesign` rejects `replicate_weights` combined with explicit `strata/psu/fpc`). diff --git a/tests/test_chaisemartin_dhaultfoeuille.py b/tests/test_chaisemartin_dhaultfoeuille.py index e09c185c..086aaa79 100644 --- a/tests/test_chaisemartin_dhaultfoeuille.py +++ b/tests/test_chaisemartin_dhaultfoeuille.py @@ -222,6 +222,48 @@ def test_convenience_function_matches_class(self): assert results_class.overall_att == pytest.approx(results_fn.overall_att) assert results_class.overall_se == pytest.approx(results_fn.overall_se) + def test_convenience_function_routes_paths_of_interest_to_init(self): + """`paths_of_interest` is an __init__ kwarg; the convenience helper + must split it out of `**kwargs` rather than letting it fall through + to fit() (which would raise TypeError). Regression for the + signature-derived split.""" + df = _by_path_three_path_data() + results_class = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)], + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + r_class = results_class.fit( + df, + outcome="outcome", + group="group", + time="period", + treatment="treatment", + L_max=3, + ) + r_fn = chaisemartin_dhaultfoeuille( + df, + outcome="outcome", + group="group", + time="period", + treatment="treatment", + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)], + twfe_diagnostic=False, + seed=42, + L_max=3, + ) + # Both surfaces produce identical per-path effects. + assert list(r_fn.path_effects.keys()) == list(r_class.path_effects.keys()) + for path in r_fn.path_effects: + for l_h, vals in r_fn.path_effects[path]["horizons"].items(): + assert vals["effect"] == pytest.approx( + r_class.path_effects[path]["horizons"][l_h]["effect"] + ) + def test_minimal_computation_path(self): # Disable everything optional; verify still works data = generate_reversible_did_data(n_groups=30, n_periods=4, seed=1) @@ -3819,31 +3861,6 @@ def test_forbids_phase3_fit_kwargs(self, fit_kwargs, msg): **fit_kwargs, ) - def test_forbids_non_binary_treatment(self): - # Continuous-dose treatment — detected inside the cell aggregator. - rng = np.random.default_rng(0) - rows = [] - for g in range(1, 7): - for t in range(4): - d = float(g % 3) if (g <= 3 and t >= 1) else 0.0 - y = d + rng.normal() - rows.append({"group": g, "period": t, "treatment": d, "outcome": y}) - data = pd.DataFrame(rows) - with warnings.catch_warnings(): - warnings.simplefilter("ignore", UserWarning) - est = ChaisemartinDHaultfoeuille( - drop_larger_lower=False, by_path=3, twfe_diagnostic=False - ) - with pytest.raises(NotImplementedError, match="non-binary"): - est.fit( - data, - outcome="outcome", - group="group", - time="period", - treatment="treatment", - L_max=2, - ) - class TestByPathBehavior: """Path enumeration, ranking, and result dict shape.""" @@ -6621,7 +6638,7 @@ def _add(group, treatment_path): str(w.message) for w in caught if issubclass(w.category, UserWarning) - and "by_path + controls" in str(w.message) + and "+ controls" in str(w.message) and "multi-baseline" not in str(w.message).lower() or ( issubclass(w.category, UserWarning) @@ -7562,7 +7579,7 @@ def _add(group, treatment_path): str(w.message) for w in caught if issubclass(w.category, UserWarning) - and "by_path + trends_linear" in str(w.message) + and "+ trends_linear" in str(w.message) and "switcher baselines" in str(w.message) ] assert deviation_msgs, ( @@ -7592,7 +7609,7 @@ def test_single_baseline_panel_does_not_emit_r_deviation_warning(self): str(w.message) for w in caught if issubclass(w.category, UserWarning) - and "by_path + trends_linear" in str(w.message) + and "+ trends_linear" in str(w.message) and "switcher baselines" in str(w.message) ] assert not deviation_msgs, ( @@ -7650,7 +7667,7 @@ def _add(group, treatment_path): str(w.message) for w in caught if issubclass(w.category, UserWarning) - and "by_path + trends_linear" in str(w.message) + and "+ trends_linear" in str(w.message) and "F_g=3" in str(w.message) ] assert boundary_msgs, ( @@ -7701,7 +7718,7 @@ def test_single_baseline_heterogeneous_F_g_does_not_warn(self): str(w.message) for w in caught if issubclass(w.category, UserWarning) - and "by_path + trends_linear" in str(w.message) + and "+ trends_linear" in str(w.message) and "switcher baselines" in str(w.message) ] assert not deviation_msgs, ( @@ -8100,3 +8117,998 @@ def test_per_path_placebos_with_trends_nonparam_bootstrap_inference(self): "trends_nonparam restriction; set_ids may not be reaching " "the per-path placebo bootstrap collector." ) + + +# --------------------------------------------------------------------------- +# Wave 3 #8: by_path + non-binary integer treatment +# --------------------------------------------------------------------------- + + +def _by_path_data_with_non_binary_treatment(seed: int = 44) -> pd.DataFrame: + """Multi-path single-baseline panel with integer-coded D in {0, 1, 2}. + + 13-period panel, 78 switchers across 3 non-binary paths + (`(0, 1, 1, 1)`, `(0, 2, 2, 2)`, `(0, 1, 2, 2)`), F_g spread + starting at 4 to keep clear of any pre-window F_g boundary cases, + plus 20 never-treated (D=0) and 20 always-treated (D=2) controls. + Outcome shifts proportionally to D so per-path effects are + distinguishable. + """ + rng = np.random.default_rng(seed) + n_periods = 13 + target_paths = [ + (0, 1, 1, 1), # path 1, low-dose sustained + (0, 2, 2, 2), # path 2, high-dose sustained + (0, 1, 2, 2), # path 3, ramp-up + ] + fg_path_counts = [ + (4, 0, 18), (5, 0, 14), # path 1 = 32 + (6, 1, 14), (7, 1, 12), # path 2 = 26 + (8, 2, 12), (9, 2, 8), # path 3 = 20 + ] + rows = [] + g_id = 0 + for F_g, path_idx, count in fg_path_counts: + target = target_paths[path_idx] + L_max = 3 + for _ in range(count): + D_row = [0] * n_periods + for j in range(L_max + 1): + D_row[F_g - 1 + j] = target[j] + for t in range(F_g + L_max, n_periods): + D_row[t] = target[L_max] + for t, d in enumerate(D_row): + rows.append({"group": g_id, "period": t, "treatment": d}) + g_id += 1 + for _ in range(20): + for t in range(n_periods): + rows.append({"group": g_id, "period": t, "treatment": 0}) + g_id += 1 + for _ in range(20): + for t in range(n_periods): + rows.append({"group": g_id, "period": t, "treatment": 2}) + g_id += 1 + df = pd.DataFrame(rows) + n_groups = df["group"].nunique() + group_fe = rng.normal(0, 2.0, size=n_groups) + df["outcome"] = ( + 10.0 + + group_fe[df["group"].values] + + 0.1 * df["period"].values + + 1.5 * df["treatment"].values + + rng.normal(0, 0.5, size=len(df)) + ) + return df + + +class TestByPathNonBinary: + """Wave 3 #8: ``by_path`` + non-binary integer treatment. + + The previous gate at ``chaisemartin_dhaultfoeuille.py:1870`` rejected + non-binary treatment + by_path. After PR Wave 3 #8+#9 it is replaced + by a D-integer validation: integer-coded D (D in Z) is supported + and continuous D (D=1.5 etc.) raises ValueError. + """ + + def test_no_longer_raises_on_non_binary(self): + """by_path=2 + D in {0,1,2} fits without raising.""" + df = _by_path_data_with_non_binary_treatment() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=2, twfe_diagnostic=False, seed=42 + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert res.path_effects is not None + assert len(res.path_effects) == 2 + + def test_non_integer_D_raises(self): + """D values containing 1.5 raise ValueError.""" + df = _by_path_data_with_non_binary_treatment() + # Cast `treatment` to float BEFORE assigning 1.5: on pandas >= 2.x + # `df.loc[mask, "treatment"] = 1.5` raises `TypeError: Invalid + # value '1.5' for dtype 'int64'` outright instead of silently + # coercing. We want to inject a continuous value to test the + # estimator's D-integer guard, not exercise pandas dtype coercion. + df["treatment"] = df["treatment"].astype(float) + mask = (df["group"] == 0) & (df["period"] == 4) + df.loc[mask, "treatment"] = 1.5 + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=2, twfe_diagnostic=False, seed=42 + ) + with pytest.raises(ValueError, match="integer-coded"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + + def test_negative_integer_D_supported(self): + """D in {-1, 0, 1} fits and produces correct path tuples.""" + rng = np.random.default_rng(45) + rows = [] + n_periods = 8 + # 30 switchers on path (0, -1, -1, -1), F_g=4 + for g in range(30): + for t in range(n_periods): + d = -1 if t >= 3 else 0 + rows.append({"group": g, "period": t, "treatment": d}) + # 30 switchers on path (0, 1, 1, 1), F_g=4 + for g in range(30, 60): + for t in range(n_periods): + d = 1 if t >= 3 else 0 + rows.append({"group": g, "period": t, "treatment": d}) + # 20 never-treated controls + for g in range(60, 80): + for t in range(n_periods): + rows.append({"group": g, "period": t, "treatment": 0}) + df = pd.DataFrame(rows) + df["outcome"] = ( + 10.0 + + df["group"].values * 0.1 + + 0.1 * df["period"].values + + 2.0 * df["treatment"].values + + rng.normal(0, 0.5, size=len(df)) + ) + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=2, twfe_diagnostic=False, seed=42 + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert res.path_effects is not None + path_keys = set(res.path_effects.keys()) + # Both paths should appear; the negative-D path must contain -1 + assert (0, -1, -1, -1) in path_keys + assert (0, 1, 1, 1) in path_keys + + def test_path_effects_present_under_non_binary(self): + """path_effects populated; tuple keys are non-binary.""" + df = _by_path_data_with_non_binary_treatment() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=3, twfe_diagnostic=False, seed=42 + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert res.path_effects is not None + assert len(res.path_effects) == 3 + # At least one path key should contain a 2 (non-binary marker). + any_non_binary = any(2 in p for p in res.path_effects.keys()) + assert any_non_binary, ( + f"Expected at least one non-binary integer in path keys, " + f"got {list(res.path_effects.keys())}" + ) + for path, entry in res.path_effects.items(): + for l_h, vals in entry["horizons"].items(): + assert np.isfinite(vals["effect"]), ( + f"path={path} l={l_h}: effect not finite" + ) + + def test_per_period_effects_unaffected_by_non_binary_by_path(self): + """per_period_effects is unchanged by the by_path lift.""" + df = _by_path_data_with_non_binary_treatment() + est_no = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=None, twfe_diagnostic=False, seed=42 + ) + est_bp = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=3, twfe_diagnostic=False, seed=42 + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res_no = est_no.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + res_bp = est_bp.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + # per_period_effects is computed before path enumeration; bit-identical. + for t in res_no.per_period_effects.keys(): + for k in ("did_plus_t", "did_minus_t"): + v_no = res_no.per_period_effects[t].get(k) + v_bp = res_bp.per_period_effects[t].get(k) + if v_no is None or not np.isfinite(v_no): + continue + assert np.isclose(v_no, v_bp, atol=1e-14, rtol=1e-14), ( + f"per_period_effects[{t}][{k}]: {v_no} != {v_bp}" + ) + + def test_to_dataframe_by_path_with_non_binary(self): + """level='by_path' DataFrame includes non-binary path-tuple labels.""" + df = _by_path_data_with_non_binary_treatment() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=3, twfe_diagnostic=False, seed=42 + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + out = res.to_dataframe(level="by_path") + assert len(out) > 0 + # path column should contain non-binary integer tuples + any_non_binary = any(2 in p for p in out["path"].unique()) + assert any_non_binary + + def test_continuous_D_without_by_path_unaffected(self): + """Continuous D + by_path=None / paths_of_interest=None: no new gate fires.""" + rng = np.random.default_rng(46) + rows = [] + # 30 switchers with continuous D in {0, 1.5} + for g in range(30): + for t in range(8): + d = 1.5 if t >= 3 else 0.0 + rows.append({"group": g, "period": t, "treatment": d}) + for g in range(30, 50): + for t in range(8): + rows.append({"group": g, "period": t, "treatment": 0.0}) + df = pd.DataFrame(rows) + df["outcome"] = ( + 10.0 + + df["group"].values * 0.1 + + 0.1 * df["period"].values + + 2.0 * df["treatment"].values + + rng.normal(0, 0.5, size=len(df)) + ) + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, twfe_diagnostic=False, seed=42 + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + # No by_path; the new D-integer validation does not fire. + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=2, + ) + assert np.isfinite(res.overall_att) + + @pytest.mark.slow + def test_bootstrap_with_non_binary_finite_se(self): + """Bootstrap SE finite on every path under non-binary D.""" + df = _by_path_data_with_non_binary_treatment() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=3, n_bootstrap=200, + twfe_diagnostic=False, seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + for path, entry in res.path_effects.items(): + for l_h, vals in entry["horizons"].items(): + assert np.isfinite(vals["se"]), ( + f"path={path} l={l_h}: bootstrap SE not finite" + ) + + @pytest.mark.slow + def test_per_path_placebos_with_non_binary_present(self): + """path_placebo_event_study populated under non-binary + placebo.""" + df = _by_path_data_with_non_binary_treatment() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=3, placebo=True, + twfe_diagnostic=False, seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert res.path_placebo_event_study is not None + assert len(res.path_placebo_event_study) > 0 + # At least one (path, -l) entry has a finite point estimate. + any_finite = False + for path, lags in res.path_placebo_event_study.items(): + for lag, vals in lags.items(): + if np.isfinite(vals["effect"]): + any_finite = True + break + if any_finite: + break + assert any_finite + + @pytest.mark.slow + def test_sup_t_bands_with_non_binary_finite_crit(self): + """Per-path sup-t crit_value finite under non-binary D.""" + df = _by_path_data_with_non_binary_treatment() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=3, n_bootstrap=400, + twfe_diagnostic=False, seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert res.path_sup_t_bands is not None + # At least one path passed the strict-majority gate. + finite_crits = [ + entry["crit_value"] + for entry in res.path_sup_t_bands.values() + if np.isfinite(entry["crit_value"]) + ] + assert len(finite_crits) > 0, ( + "Expected at least one finite crit_value under non-binary D + bootstrap" + ) + + +# --------------------------------------------------------------------------- +# Wave 3 #9: paths_of_interest (Python-only API, mutex with by_path) +# --------------------------------------------------------------------------- + + +class TestPathsOfInterest: + """``paths_of_interest`` user-specified path subset. + + Validation, mutex with ``by_path``, behavior, and cross-feature + composition. No R parity (R has no list-based path selection). + """ + + # ---- __init__ validation ---- + + def test_invalid_type_raises(self): + with pytest.raises(ValueError, match="paths_of_interest must be"): + ChaisemartinDHaultfoeuille(paths_of_interest="not a list") + + def test_empty_list_raises(self): + with pytest.raises(ValueError, match="must be non-empty"): + ChaisemartinDHaultfoeuille(paths_of_interest=[]) + + def test_non_tuple_path_raises(self): + with pytest.raises(ValueError, match=r"paths_of_interest\[0\]"): + ChaisemartinDHaultfoeuille(paths_of_interest=[{0, 1, 1, 1}]) + + def test_non_int_element_raises(self): + with pytest.raises(ValueError, match="must be an int"): + ChaisemartinDHaultfoeuille( + paths_of_interest=[(0, "a", 1, 1)] + ) + + def test_bool_element_raises(self): + with pytest.raises(ValueError, match="must be an int"): + ChaisemartinDHaultfoeuille( + paths_of_interest=[(False, True, True, True)] + ) + + def test_np_bool_element_raises(self): + with pytest.raises(ValueError, match="must be an int"): + ChaisemartinDHaultfoeuille( + paths_of_interest=[(np.bool_(True), 0, 0, 0)] + ) + + def test_np_integer_accepted_canonicalized(self): + est = ChaisemartinDHaultfoeuille( + paths_of_interest=[(np.int64(0), np.int32(1), 1, 1)] + ) + # Canonicalized to Python int tuples. + assert est.paths_of_interest == [(0, 1, 1, 1)] + for v in est.paths_of_interest[0]: + assert type(v) is int + + def test_mixed_lengths_raise(self): + with pytest.raises(ValueError, match="mixed lengths"): + ChaisemartinDHaultfoeuille( + paths_of_interest=[(0, 1), (0, 1, 1, 1)] + ) + + def test_mutex_with_by_path_raises(self): + with pytest.raises(ValueError, match="mutually exclusive"): + ChaisemartinDHaultfoeuille( + by_path=2, paths_of_interest=[(0, 1, 1, 1)] + ) + + def test_set_params_re_validates_mutex_by_path_added(self): + est = ChaisemartinDHaultfoeuille(paths_of_interest=[(0, 1, 1, 1)]) + with pytest.raises(ValueError, match="mutually exclusive"): + est.set_params(by_path=2) + + def test_set_params_re_validates_mutex_poi_added(self): + """Reciprocal: construct with by_path, set_params(paths_of_interest=...).""" + est = ChaisemartinDHaultfoeuille(by_path=2) + with pytest.raises(ValueError, match="mutually exclusive"): + est.set_params(paths_of_interest=[(0, 1, 1, 1)]) + + def test_set_params_failed_validation_is_transactional(self): + """A failed `set_params()` must leave estimator state unchanged + (regression for R5 P2 finding: prior implementation mutated + before validation, leaving both selectors populated when the + mutex check raised, which `fit()` then silently consumed).""" + est = ChaisemartinDHaultfoeuille(paths_of_interest=[(0, 1, 1, 1)]) + # Capture pre-call state. + before = est.get_params() + # Mutex violation: by_path AND paths_of_interest both non-None. + with pytest.raises(ValueError, match="mutually exclusive"): + est.set_params(by_path=2) + # Post-failure state is rolled back to pre-call. + after = est.get_params() + assert after == before, ( + f"set_params() rollback failed: by_path={after['by_path']}, " + f"paths_of_interest={after['paths_of_interest']}" + ) + # Subsequent valid set_params() succeeds against rolled-back state. + est.set_params(by_path=2, paths_of_interest=None) + params = est.get_params() + assert params["by_path"] == 2 + assert params["paths_of_interest"] is None + + def test_get_params_includes_paths_of_interest(self): + est = ChaisemartinDHaultfoeuille( + paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)] + ) + params = est.get_params() + assert "paths_of_interest" in params + assert params["paths_of_interest"] == [(0, 1, 1, 1), (0, 1, 0, 0)] + + def test_canonicalized_duplicates_dedup_warn(self): + """Cross-numeric-type duplicates collapse and warn at fit-time.""" + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[ + (np.int64(0), 1, 1, 1), + (0, 1, 1, 1), + ], + twfe_diagnostic=False, + seed=42, + ) + # Canonicalization at __init__ produces two identical Python int + # tuples; the dedup warning fires inside _enumerate_treatment_paths. + with pytest.warns(UserWarning, match="duplicate path"): + with warnings.catch_warnings(): + warnings.simplefilter("default", UserWarning) + est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + + # ---- fit-time validation ---- + + def test_wrong_length_raises_at_fit(self): + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1)], # length 3, expected L_max+1=4 + twfe_diagnostic=False, + seed=42, + ) + with pytest.raises(ValueError, match="length L_max"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + + def test_paths_of_interest_requires_L_max(self): + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1)], + twfe_diagnostic=False, + seed=42, + ) + with pytest.raises(ValueError, match="requires L_max"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", + ) + + def test_paths_of_interest_requires_drop_larger_lower_false(self): + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + paths_of_interest=[(0, 1, 1, 1)], + twfe_diagnostic=False, + seed=42, + ) + with pytest.raises(ValueError, match="drop_larger_lower=False"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + + # ---- behavior ---- + + def test_paths_of_interest_selects_user_paths(self): + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)], + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert set(res.path_effects.keys()) == {(0, 1, 1, 1), (0, 1, 0, 0)} + + def test_paths_of_interest_preserves_user_order(self): + df = _by_path_three_path_data() + # Order intentionally NOT frequency-ranked: lower-rank path first. + user_order = [(0, 1, 0, 0), (0, 1, 1, 1)] + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=user_order, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + # Insertion order preserved. + assert list(res.path_effects.keys()) == user_order + + def test_paths_of_interest_order_preserved_in_to_dataframe(self): + """`to_dataframe(level="by_path")` must iterate in insertion + order so user-specified `paths_of_interest` order is preserved + across reporting surfaces (regression for R1 P3 + maintainability finding).""" + df = _by_path_three_path_data() + # Order intentionally NOT frequency-ranked: lower-rank path first. + user_order = [(0, 1, 0, 0), (0, 1, 1, 1)] + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=user_order, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + out = res.to_dataframe(level="by_path") + # First-occurrence order of the path column matches user order. + first_seen = [] + seen = set() + for p in out["path"]: + if p not in seen: + seen.add(p) + first_seen.append(p) + assert first_seen == user_order + + def test_paths_of_interest_order_preserved_in_summary(self): + """`summary()` must render paths in insertion order so + user-specified `paths_of_interest` order is preserved + (regression for R1 P3 maintainability finding).""" + df = _by_path_three_path_data() + user_order = [(0, 1, 0, 0), (0, 1, 1, 1)] + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=user_order, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + text = res.summary() + # Find the ordering of the user paths as they appear in summary. + idx_first = text.index("Path (0, 1, 0, 0)") + idx_second = text.index("Path (0, 1, 1, 1)") + assert idx_first < idx_second, ( + f"Summary did not preserve user-specified order. " + f"`(0, 1, 0, 0)` at idx={idx_first}, " + f"`(0, 1, 1, 1)` at idx={idx_second}" + ) + + def test_paths_of_interest_frequency_rank_is_true_frequency(self): + """`frequency_rank` must reflect descending count, NOT user-list + order. Regression for the R0 P2 finding: previously the rank + field was assigned from `enumerate(selected_paths)` which gave + user-selection order under `paths_of_interest`.""" + df = _by_path_three_path_data() + # _by_path_three_path_data: (0,1,1,1) has 3 groups, (0,1,0,0) has 2, + # (0,1,1,0) has 1. User passes the lowest-frequency path first. + user_order = [(0, 1, 0, 0), (0, 1, 1, 1)] + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=user_order, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + # (0,1,1,1) has higher frequency → rank 1 + # (0,1,0,0) has lower frequency → rank 2 + assert res.path_effects[(0, 1, 1, 1)]["frequency_rank"] == 1 + assert res.path_effects[(0, 1, 0, 0)]["frequency_rank"] == 2 + + def test_paths_of_interest_all_unobserved_summary_distinct_text(self): + """When every path in `paths_of_interest` is unobserved, the + empty-state `summary()` block must render the + paths_of_interest-specific text rather than the generic + "no observed paths have a complete window" text (regression + for R4 P3 maintainability finding).""" + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(1, 1, 1, 1), (1, 0, 1, 0)], # both unobserved + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + text = res.summary() + assert "Every path in paths_of_interest was unobserved" in text, ( + f"summary() did not render the paths_of_interest-specific " + f"empty-state text. Got:\n{text}" + ) + # And the generic by_path-only wording must NOT appear in this case. + assert "No observed paths have a complete" not in text + + def test_paths_of_interest_all_unobserved_emits_distinct_warning(self): + """When every path in `paths_of_interest` is unobserved, + the empty-state warning should mention `paths_of_interest` + explicitly rather than the generic `by_path={n}` text + (regression for R2 P3 maintainability finding).""" + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(1, 1, 1, 1), (1, 0, 1, 0)], # both unobserved + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(record=True) as recorded: + warnings.simplefilter("always", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + # The summary empty-state warning mentions paths_of_interest, not by_path + empty_state_warnings = [ + str(w.message) + for w in recorded + if "paths_of_interest was requested but every" in str(w.message) + ] + assert len(empty_state_warnings) >= 1, ( + f"Expected paths_of_interest-specific empty-state warning, " + f"got: {[str(w.message) for w in recorded]}" + ) + # And the result is empty dict, not None + assert res.path_effects == {} + + def test_unobserved_path_warns_and_omits(self): + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[ + (0, 1, 1, 1), + (1, 1, 1, 1), # not observed (no group has D=1 at F_g-1) + ], + twfe_diagnostic=False, + seed=42, + ) + with pytest.warns(UserWarning, match="zero observed"): + with warnings.catch_warnings(): + warnings.simplefilter("default", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert (1, 1, 1, 1) not in res.path_effects + assert (0, 1, 1, 1) in res.path_effects + + @pytest.mark.slow + def test_paths_of_interest_with_non_binary_D(self): + df = _by_path_data_with_non_binary_treatment() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 2, 2, 2)], + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert set(res.path_effects.keys()) == {(0, 2, 2, 2)} + for l_h, vals in res.path_effects[(0, 2, 2, 2)]["horizons"].items(): + assert np.isfinite(vals["effect"]) + + # ---- cross-feature composition (review MEDIUM #3) ---- + + @pytest.mark.slow + def test_paths_of_interest_with_controls(self): + df = _by_path_data_with_non_binary_treatment().copy() + # Inject a single-baseline control so multi-baseline warning doesn't fire + df["X1"] = np.random.default_rng(0).normal(size=len(df)) + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1), (0, 2, 2, 2)], + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, controls=["X1"], + ) + assert len(res.path_effects) == 2 + for path, entry in res.path_effects.items(): + for l_h, vals in entry["horizons"].items(): + assert np.isfinite(vals["effect"]), f"path={path} l={l_h}" + + @pytest.mark.slow + def test_paths_of_interest_with_trends_linear(self): + df = _by_path_data_with_trends_linear() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)], + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, trends_linear=True, + ) + assert len(res.path_effects) == 2 + # path_cumulated_event_study populated under trends_linear. + assert res.path_cumulated_event_study is not None + assert set(res.path_cumulated_event_study.keys()) == set(res.path_effects.keys()) + + @pytest.mark.slow + def test_paths_of_interest_with_trends_nonparam(self): + df = _by_path_data_with_trends_nonparam() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)], + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + trends_nonparam="state", + ) + assert len(res.path_effects) == 2 + for path, entry in res.path_effects.items(): + for l_h, vals in entry["horizons"].items(): + assert np.isfinite(vals["effect"]), f"path={path} l={l_h}" + + @pytest.mark.slow + def test_paths_of_interest_non_binary_bootstrap_placebo(self): + """Quadruple-combo: paths_of_interest + non-binary D + bootstrap + + placebo. Asserts (i) selector restricts path_effects to the + requested paths, (ii) bootstrap SE finite on every (path, + horizon) for the selected paths, (iii) per-path placebo + populated, (iv) at least one selected path has a finite sup-t + crit_value and the corresponding `cband_conf_int` is non-NaN.""" + df = _by_path_data_with_non_binary_treatment() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 2, 2, 2), (0, 1, 1, 1)], + n_bootstrap=400, + placebo=True, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert set(res.path_effects.keys()) == {(0, 2, 2, 2), (0, 1, 1, 1)} + # 1. analytical / bootstrap on path_effects (SE finite) + for path in res.path_effects: + for l_h, vals in res.path_effects[path]["horizons"].items(): + assert np.isfinite(vals["effect"]), f"path={path} l={l_h}" + assert np.isfinite(vals["se"]), f"path={path} l={l_h}" + # 2. per-path placebo + assert res.path_placebo_event_study is not None + assert (0, 2, 2, 2) in res.path_placebo_event_study + assert (0, 1, 1, 1) in res.path_placebo_event_study + # 3. per-path sup-t bands: at least one selected path passes the + # strict-majority gate with a finite crit_value AND the + # corresponding cband_conf_int entries on path_effects are + # non-NaN tuples (vacuous "is not None" check rejected by R6). + assert res.path_sup_t_bands is not None + finite_crit_paths = [ + p + for p, entry in res.path_sup_t_bands.items() + if np.isfinite(entry.get("crit_value", np.nan)) + ] + assert len(finite_crit_paths) >= 1, ( + f"Expected >=1 selected path with finite sup-t crit; " + f"got path_sup_t_bands={res.path_sup_t_bands}" + ) + # cband_conf_int populated for positive horizons of finite-crit paths. + for path in finite_crit_paths: + for l_h in range(1, 4): + cband = res.path_effects[path]["horizons"][l_h].get( + "cband_conf_int" + ) + assert cband is not None, f"path={path} l={l_h}: cband missing" + lo, hi = cband + assert np.isfinite(lo) and np.isfinite(hi), ( + f"path={path} l={l_h}: cband endpoints not finite " + f"(lo={lo}, hi={hi})" + ) + + @pytest.mark.slow + def test_paths_of_interest_trends_linear_bootstrap_placebo(self): + """`paths_of_interest + trends_linear=True + n_bootstrap > 0 + + placebo=True`: assert (i) selector restricts the path set, (ii) + per-path bootstrap SE on `path_effects` finite, (iii) + post-bootstrap `path_cumulated_event_study` populated for the + same paths (mirrors global `linear_trends_effects` cumulation), + (iv) per-path placebo populated. Regression for the R6 P1 + Cartesian-product gap.""" + df = _by_path_data_with_trends_linear() + user_paths = [(0, 1, 1, 1), (0, 1, 1, 0)] + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=user_paths, + n_bootstrap=200, + placebo=True, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, trends_linear=True, + ) + # (i) selector restricts the path set + assert set(res.path_effects.keys()) == set(user_paths) + # (ii) per-path bootstrap SE finite on event study + for path in res.path_effects: + for l_h, vals in res.path_effects[path]["horizons"].items(): + assert np.isfinite(vals["se"]), ( + f"path={path} l={l_h}: bootstrap SE not finite" + ) + # (iii) post-bootstrap path_cumulated_event_study populated for + # the same paths AND derived from the post-bootstrap per-horizon + # SEs (cumulated SE = sum of post-bootstrap component SEs). + assert res.path_cumulated_event_study is not None + assert set(res.path_cumulated_event_study.keys()) == set(user_paths) + for path in user_paths: + assert len(res.path_cumulated_event_study[path]) > 0 + for l_h, vals in res.path_cumulated_event_study[path].items(): + assert np.isfinite(vals["effect"]), ( + f"path={path} l={l_h}: cumulated effect not finite" + ) + assert np.isfinite(vals["se"]), ( + f"path={path} l={l_h}: cumulated SE not finite" + ) + # (iv) per-path placebo populated + assert res.path_placebo_event_study is not None + assert set(res.path_placebo_event_study.keys()) == set(user_paths) + + @pytest.mark.slow + def test_paths_of_interest_trends_nonparam_bootstrap_placebo(self): + """`paths_of_interest + trends_nonparam="state" + placebo=True + + n_bootstrap > 0`: assert the selector + set_ids flow through + the four per-path collectors (`_compute_path_effects`, + `_compute_path_placebos`, `_collect_path_bootstrap_inputs`, + `_collect_path_placebo_bootstrap_inputs`) and the resulting + public surfaces are populated. Regression for the R6 P1 + Cartesian-product gap.""" + df = _by_path_data_with_trends_nonparam() + user_paths = [(0, 1, 1, 1), (0, 1, 1, 0)] + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=user_paths, + n_bootstrap=200, + placebo=True, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + trends_nonparam="state", + ) + # Selector restriction + assert set(res.path_effects.keys()) == set(user_paths) + # Per-path bootstrap SE finite on event study (bootstrap + # collector + set_ids reached path_effects) + for path in res.path_effects: + for l_h, vals in res.path_effects[path]["horizons"].items(): + assert np.isfinite(vals["effect"]), f"path={path} l={l_h}" + assert np.isfinite(vals["se"]), f"path={path} l={l_h}" + # Per-path placebo populated and bootstrap SE finite + # (placebo collector + set_ids reached path_placebo_event_study) + assert res.path_placebo_event_study is not None + assert set(res.path_placebo_event_study.keys()) == set(user_paths) + any_finite_placebo_se = False + for path in user_paths: + for lag, vals in res.path_placebo_event_study[path].items(): + if np.isfinite(vals["effect"]) and np.isfinite(vals["se"]): + any_finite_placebo_se = True + break + if any_finite_placebo_se: + break + assert any_finite_placebo_se, ( + "No finite (effect, SE) pair on per-path placebo surface " + "under paths_of_interest + trends_nonparam + bootstrap" + ) + + # ---- single-surface inheritance (slow) ---- + + @pytest.mark.slow + def test_bootstrap_with_paths_of_interest_finite_se(self): + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)], + n_bootstrap=200, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + for path, entry in res.path_effects.items(): + for l_h, vals in entry["horizons"].items(): + assert np.isfinite(vals["se"]), f"path={path} l={l_h}: bootstrap SE not finite" + + @pytest.mark.slow + def test_per_path_placebos_with_paths_of_interest_present(self): + df = _by_path_three_path_data() + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, + paths_of_interest=[(0, 1, 1, 1), (0, 1, 0, 0)], + placebo=True, + twfe_diagnostic=False, + seed=42, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + res = est.fit( + df, outcome="outcome", group="group", time="period", + treatment="treatment", L_max=3, + ) + assert res.path_placebo_event_study is not None + assert len(res.path_placebo_event_study) == 2 diff --git a/tests/test_chaisemartin_dhaultfoeuille_parity.py b/tests/test_chaisemartin_dhaultfoeuille_parity.py index 0a916cc3..ef1fbe4b 100644 --- a/tests/test_chaisemartin_dhaultfoeuille_parity.py +++ b/tests/test_chaisemartin_dhaultfoeuille_parity.py @@ -1186,3 +1186,136 @@ def test_parity_multi_path_reversible_by_path_trends_nonparam( f"path={path_key} h={h} SE: " f"py={py_se:.4f} vs r={r_se:.4f}" ) + + +class TestDCDHDynRParityByPathNonBinary: + """ + Parity tests for ``by_path`` + non-binary integer treatment against + R DIDmultiplegtDYN 2.3.3. + + Wave 3 #8 lift. R's ``did_multiplegt_dyn(..., by_path=k)`` accepts + discrete non-binary treatment naturally; the per-path string + encoding (``"0,2,2,2"``) maps bit-for-bit to Python's tuple key + ``(0, 2, 2, 2)`` for D in {0..9}. R's + ``substr(path_index$path, 1, 1)`` baseline-derivation (in + ``did_multiplegt_by_path``) breaks for D >= 10 (it captures only + the first character of the comma-separated path string), but + Python's tuple-key matching is correct in that regime — the + parity scenario stays in {0, 1, 2} so the R bug does not + interfere. + + The fixture mirrors Scenario 17's single-baseline custom DGP: + 78 switchers across 3 non-binary paths, 20 never-treated and 20 + always-treated controls, F_g >= 4. Per-path point estimates + match R bit-exactly. Per-path SE inherits the documented + cross-path cohort-sharing deviation (Phase 2 envelope, ~5% rtol + on this scenario; SE_RTOL=0.15 for headroom). + """ + + POINT_RTOL = 1e-9 + # Near-zero placebo (e.g. -0.003) amplifies relative difference; + # use atol=1e-9 as a fallback so the cumulated relative+absolute + # `pytest.approx` envelope matches bit-exact event horizons and + # small-value placebo points consistently. + POINT_ATOL = 1e-9 + SE_RTOL = 0.15 + + def _path_key_from_r_label(self, r_label: str): + return tuple(int(x) for x in r_label.split(",")) + + def test_parity_multi_path_reversible_by_path_non_binary( + self, golden_values + ): + """3-path case with D in {0, 1, 2}: by_path=3, placebo=1.""" + import math + import warnings + + scenario = golden_values.get( + "multi_path_reversible_by_path_non_binary" + ) + if scenario is None: + pytest.skip( + "scenario 'multi_path_reversible_by_path_non_binary' " + "not in golden values" + ) + + df = _golden_to_df(scenario["data"]) + est = ChaisemartinDHaultfoeuille( + drop_larger_lower=False, by_path=3, placebo=True + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + results = est.fit( + df, + outcome="outcome", + group="group", + time="period", + treatment="treatment", + L_max=3, + ) + + r_by_path = scenario["results"]["by_path"] + assert results.path_effects is not None + + py_keys = set(results.path_effects.keys()) + r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path} + assert py_keys == r_keys, ( + f"Path-set mismatch.\n" + f" Python only: {py_keys - r_keys}\n" + f" R only: {r_keys - py_keys}" + ) + + for r_path_entry in r_by_path: + path_key = self._path_key_from_r_label(r_path_entry["path"]) + py_path = results.path_effects[path_key] + py_placebos = ( + results.path_placebo_event_study.get(path_key, {}) + if results.path_placebo_event_study is not None + else {} + ) + + assert py_path["frequency_rank"] == r_path_entry["frequency_rank"], ( + f"path={path_key}: frequency_rank mismatch " + f"py={py_path['frequency_rank']} vs r={r_path_entry['frequency_rank']}" + ) + + for h_str, r_h in r_path_entry["horizons"].items(): + h = int(h_str) + if h > 0: + assert h in py_path["horizons"], ( + f"path={path_key}: horizon {h} missing from " + f"Python path_effects" + ) + py_h = py_path["horizons"][h] + else: + assert h in py_placebos, ( + f"path={path_key}: placebo {h} missing from " + f"Python path_placebo_event_study" + ) + py_h = py_placebos[h] + + assert py_h["n_obs"] == int(r_h["n_switchers"]), ( + f"path={path_key} h={h}: switcher-count mismatch " + f"py={py_h['n_obs']} vs r={int(r_h['n_switchers'])}" + ) + + assert py_h["effect"] == pytest.approx( + r_h["effect"], rel=self.POINT_RTOL, abs=self.POINT_ATOL + ), ( + f"path={path_key} h={h}: " + f"py={py_h['effect']:.6f} vs r={r_h['effect']:.6f}" + ) + + py_se = py_h["se"] + r_se = r_h["se"] + py_finite_positive = math.isfinite(py_se) and py_se > 0.0 + r_finite_positive = math.isfinite(r_se) and r_se > 0.0 + assert py_finite_positive == r_finite_positive, ( + f"path={path_key} h={h} SE state mismatch " + f"(py_se={py_se}, r_se={r_se})" + ) + if py_finite_positive and r_finite_positive: + assert py_se == pytest.approx(r_se, rel=self.SE_RTOL), ( + f"path={path_key} h={h} SE: " + f"py={py_se:.4f} vs r={r_se:.4f}" + )