diff --git a/CHANGELOG.md b/CHANGELOG.md index ed8ab360..c03efe6a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **`ContinuousDiD` lowest-dose-as-control** (`control_group="lowest_dose"`, CGBS 2024 Remark 3.1) for + settings with no untreated group (`P(D=0) = 0`): the lowest-dose group `d_L` becomes the comparison + and the estimand is `ATT(d) − ATT(d_L)` (with `ATT(d_L) = 0` the omitted reference). It is a + control-group swap — the entire linear influence-function / bootstrap / event-study / survey + machinery is reused unchanged (`ee_control` already carries the reference-group variance) — so both + the discrete (`treatment_type="discrete"`) and continuous (B-spline) paths are supported; the + discrete ACRT backward-difference reference shifts from `0` to `d_L` + (`ACRT(d_1) = ATT(d_1)/(d_1 − d_L)`). The continuous path requires a genuine mass point at the + minimum dose (`P(D=d_L) > 0`). Results gain a `reference_dose` field (`= d_L`). Fail-closed: + never-treated units present, a singleton `d_L`, no treated dose above `d_L`, `dvals ≤ d_L`, or a + survey/subpopulation design leaving `< 2` positive-weight `d_L` units all raise; multi-cohort and + `covariates=` × `lowest_dose` raise `NotImplementedError` (deferred). The default `never_treated` / + `not_yet_treated` paths are unchanged. - **`ContinuousDiD` discrete-treatment saturated regression** (`treatment_type="discrete"`) for multi-valued / discrete dose (CGBS 2024 Eq. 4.1). Each distinct dose level gets its own effect coefficient — `ATT(d_j) = mean_{D=d_j}(ΔY) − control` (a per-level 2×2 DiD) — instead of a B-spline diff --git a/TODO.md b/TODO.md index 09ef8439..0825f762 100644 --- a/TODO.md +++ b/TODO.md @@ -29,7 +29,7 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| | `SyntheticControl` conformal (CWZ 2021) extensions: (a) one-sided / signed-`t` variants (§7); (b) covariates in the conformal proxy (`X_jt`, eqs 4/6 — current proxy is outcomes-only); (c) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies. The joint test, pointwise CIs, and average-effect CI have landed. | `conformal.py`, `synthetic_control_results.py` | CWZ-2021 | Heavy | Low | -| `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; (b) discrete-treatment saturated regression (`treatment_type="discrete"`) — **DONE**; remaining: (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0`. Also deferred: `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate — documented `NotImplementedError`); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); and multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it — currently `NotImplementedError`; single-cohort / 2-period / shared-support multi-cohort are supported). | `continuous_did.py` | CGBS-2024 | Heavy | Low | +| `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; (b) discrete-treatment saturated regression (`treatment_type="discrete"`) — **DONE**; (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0` (`control_group="lowest_dose"`) — **DONE** (discrete + continuous mass-point, single-cohort; estimand `ATT(d)−ATT(d_L)`; see REGISTRY Note #7). Remaining (all deferred `NotImplementedError`, documented): `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it); **multi-cohort `lowest_dose`** (within-cohort `d_L` reference + support-aware cross-cohort aggregation); and **`covariates=` × `lowest_dose`** (conditional-PT-relative-to-`d_L` estimand). Single-cohort / 2-period / shared-support multi-cohort are supported. | `continuous_did.py` | CGBS-2024 | Heavy | Low | | `ImputationDiD` LOO conservative-variance refinement (BJS 2024 Supp. Appendix A.9) — a finite-sample improvement to the auxiliary-model residuals reducing overfit of `tau_tilde_g` to `epsilon`. Asymptotic Theorem-3 variance is implemented and matches R `didimputation` (which also omits LOO by default). | `imputation.py` | imputation-validation | Mid | Low | | `TwoWayFixedEffects(vcov_type in {hc2, hc2_bm})` with replicate-weight designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate, which doesn't compose with the full-dummy HC2/HC2-BM build — a correct impl needs per-replicate full-dummy refit. Workaround: `hc1` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Heavy | Low | | `CallawaySantAnna` DR no-covariate per-cell SE still uses the ddof=1 plug-in `sqrt(var_t/n_t + var_c/n_c)` rather than the IF-based `sqrt(sum(phi^2))` the reg/ipw paths switched to (O(1/n) from R's `drdid_panel`; aggregated SEs unaffected — they already consume the IF). Switch for method-uniform per-cell/aggregated consistency. | `staggered.py::_doubly_robust` (no-cov branch) | CS-scaling | Quick | Low | diff --git a/diff_diff/continuous_did.py b/diff_diff/continuous_did.py index badbb3ca..ead13009 100644 --- a/diff_diff/continuous_did.py +++ b/diff_diff/continuous_did.py @@ -61,7 +61,13 @@ class ContinuousDiD: dvals : array-like, optional Custom dose evaluation grid. If None, uses quantile-based default. control_group : str, default="never_treated" - ``"never_treated"`` or ``"not_yet_treated"``. + ``"never_treated"``, ``"not_yet_treated"``, or ``"lowest_dose"``. + ``"lowest_dose"`` implements Remark 3.1 (CGBS 2024) for settings with no + never-treated / zero-dose units (``P(D=0) = 0``): the lowest-dose group + ``d_L`` becomes the comparison and the estimand is ``ATT(d) − ATT(d_L)``. + Requires a genuine lowest-dose group (``>= 2`` units at ``d_L``, i.e. + ``P(D=d_L) > 0``) and no never-treated units present. Single-cohort only + (multi-cohort and ``covariates=`` raise ``NotImplementedError``). anticipation : int, default=0 Number of periods of treatment anticipation. base_period : str, default="varying" @@ -131,7 +137,7 @@ class ContinuousDiD: >>> results.overall_att # doctest: +SKIP """ - _VALID_CONTROL_GROUPS = {"never_treated", "not_yet_treated"} + _VALID_CONTROL_GROUPS = {"never_treated", "not_yet_treated", "lowest_dose"} _VALID_BASE_PERIODS = {"varying", "universal"} _VALID_ESTIMATION_METHODS = {"reg", "dr", "ipw"} _VALID_TREATMENT_TYPES = {"continuous", "discrete"} @@ -210,6 +216,17 @@ def _validate_constrained_params(self) -> None: f"Invalid treatment_type: '{self.treatment_type}'. " f"Must be one of {self._VALID_TREATMENT_TYPES}." ) + if self.control_group == "lowest_dose" and self.covariates is not None: + # The covariate estimand under lowest-dose-as-control shifts to + # conditional PT *relative to d_L* (E[ΔY(0)|D=d,X] = E[ΔY(0)|d_L,X]). + # Deferred (see TODO) rather than silently estimated with the wrong + # identifying assumption. + raise NotImplementedError( + "control_group='lowest_dose' does not yet compose with covariates= " + "(the conditional-parallel-trends estimand relative to the lowest " + "dose d_L is deferred). Use covariates=None for the unconditional " + "lowest-dose fit." + ) def get_params(self) -> Dict[str, Any]: """Return estimator parameters as a dictionary.""" @@ -509,11 +526,80 @@ def fit( if self.control_group == "not_yet_treated" and n_control == 0: raise ValueError( "No never-treated (D=0) units found. With control_group='not_yet_treated', " - "dose-response curve identification requires P(D=0) > 0 " - "(Remark 3.1 in Callaway et al. is not yet implemented). " - "Add never-treated units or use a dataset with D=0 observations." + "dose-response curve identification requires P(D=0) > 0. For settings " + "with no untreated group, use control_group='lowest_dose' (Remark 3.1: " + "the lowest-dose group becomes the comparison, estimand ATT(d)-ATT(d_L)). " + "Otherwise add never-treated units or use a dataset with D=0 observations." ) + # Remark 3.1 (control_group="lowest_dose"): the lowest-dose group d_L is + # the comparison. Compute d_L ONCE here (from the treated unit doses, + # before precompute) and thread it via precomp -> every d_L-referencing + # consumer runs after this, and fit() stays config-idempotent (no fitted + # self attr). The d_L cluster (|dose - d_L| <= SATURATED_TOL) is the + # single source of truth for both the mask and the modelled dose set. + lowest_dose: Optional[float] = None + if self.control_group == "lowest_dose": + if n_control > 0: + raise ValueError( + "control_group='lowest_dose' is for settings with no never-treated " + f"units (Remark 3.1, P(D=0)=0), but {n_control} never-treated unit(s) " + "were found; they would be silently dropped. Use " + "control_group='never_treated' or 'not_yet_treated', or remove the " + "never-treated units." + ) + if len(treatment_groups) > 1: + # NOTE (deferred multi-cohort follow-up): a future multi-cohort + # lowest_dose must use a WITHIN-cohort d_L reference and a + # support-aware cross-cohort aggregation, and must exclude the d_L + # controls from the survey group/bin mass sums (which key off + # unit_cohorts==g and would otherwise double-count them). Harmless + # today because this path is fenced off here. + raise NotImplementedError( + "control_group='lowest_dose' with multiple treatment cohorts is not " + f"yet implemented ({len(treatment_groups)} cohorts found). Remark 3.1 " + "is defined for a single treatment date; use a single-cohort panel " + "(multi-period single-cohort is supported)." + ) + dose_arr = treated_unit_doses.to_numpy(dtype=float) + d_L = float(np.min(dose_arr)) + n_dL = int(np.sum(np.abs(dose_arr - d_L) <= SATURATED_TOL)) + if n_dL < 2: + msg = ( + f"control_group='lowest_dose' requires a lowest-dose *group* — a mass " + f"point at the minimum dose d_L={d_L:g} with >= 2 units (P(D=d_L) > 0), " + f"but only {n_dL} unit is at d_L. The reference group must have enough " + "units to form its own control variance; a singleton minimum is not a " + "lowest-dose group." + ) + if self.treatment_type != "discrete": + msg += ( + " On a truly continuous dose without a mass point at the minimum, " + "Remark 3.1 does not apply." + ) + raise ValueError(msg) + above = dose_arr[dose_arr - d_L > SATURATED_TOL] + if len(saturated_dose_levels(above)) < 1: + raise ValueError( + f"control_group='lowest_dose': no treated dose above the lowest dose " + f"d_L={d_L:g}. The estimand ATT(d)-ATT(d_L) needs at least one dose " + "level above d_L, but all treated units share the same dose." + ) + # A lowest modelled dose d_1 very close to d_L makes the boundary + # ACRT(d_1)=ATT(d_1)/(d_1-d_L) and its SE explode; warn (not an error). + d_1 = float(np.min(above)) + dose_span = float(np.max(dose_arr) - d_L) + if dose_span > 0 and (d_1 - d_L) < 0.01 * dose_span: + warnings.warn( + f"control_group='lowest_dose': the lowest modelled dose d_1={d_1:g} is " + f"very close to the reference d_L={d_L:g} (gap {d_1 - d_L:g}, " + f"{100 * (d_1 - d_L) / dose_span:.2g}% of the dose range); the boundary " + "ACRT(d_1)=ATT(d_1)/(d_1-d_L) and its standard error may be very large.", + UserWarning, + stacklevel=2, + ) + lowest_dose = d_L + # Re-resolve survey design on filtered df if rows were dropped # (survey arrays must align with df, not the original data) if resolved_survey is not None and len(df) < len(data): @@ -533,13 +619,53 @@ def fit( survey_weights=survey_weights, covariates=self.covariates, ) + # Thread the lowest-dose reference d_L (Remark 3.1) to the per-cell + # dose-response so it swaps the control group and shifts the discrete + # ACRT reference. None on the never/not-yet-treated paths. + precomp["lowest_dose"] = lowest_dose # Compute dvals (evaluation grid); for discrete treatment, also the - # saturated dose levels (the global basis support). + # saturated dose levels (the global basis support). Under lowest_dose the + # lowest-dose group d_L is the reference (not modelled): the basis / grid + # / levels span only the *modelled* doses strictly above d_L. all_treated_doses = precomp["dose_vector"][precomp["dose_vector"] > 0] + if lowest_dose is not None: + modelled_doses = all_treated_doses[all_treated_doses - lowest_dose > SATURATED_TOL] + if self.dvals is not None: + bad = self.dvals[self.dvals <= lowest_dose + SATURATED_TOL] + if bad.size: + raise ValueError( + f"control_group='lowest_dose': dvals contain {int(bad.size)} " + f"value(s) <= the reference dose d_L={lowest_dose:g}. d_L is the " + "omitted reference (ATT(d_L)=0 by construction); evaluate the " + "dose-response only at doses strictly above d_L." + ) + # Survey subpopulation weights could reduce the d_L control group to + # zero or a single positive-weight unit, leaving no identified + # reference to difference against: with one positive-weight d_L unit, + # mu_0 equals that unit's dY so its ee_control = w*(dY - mu_0) = 0 and + # the reference contributes zero variance (understated SE). The raw + # >= 2 guard runs before survey weights, so enforce the same effective + # >= 2 positive-weight requirement here. Fail closed. + usw0 = precomp.get("unit_survey_weights") + if usw0 is not None: + dv0 = precomp["dose_vector"] + dL_w = usw0[np.abs(dv0 - lowest_dose) <= SATURATED_TOL] + n_pos_dL = int(np.count_nonzero(dL_w > 0)) + if n_pos_dL < 2: + raise ValueError( + f"control_group='lowest_dose': the lowest-dose group d_L=" + f"{lowest_dose:g} has {n_pos_dL} positive-weight unit(s) after " + "survey/subpopulation weighting (< 2 needed for an identified " + "reference variance; a single positive-weight reference unit " + "contributes zero control-side variance). Widen the subpopulation " + "or use a different dose grid." + ) + else: + modelled_doses = all_treated_doses levels: Optional[np.ndarray] = None if self.treatment_type == "discrete": - levels = saturated_dose_levels(all_treated_doses) + levels = saturated_dose_levels(modelled_doses) if self.dvals is not None: # A saturated model can only be evaluated at observed dose # levels; reject an off-support request (no silent snapping). @@ -602,12 +728,13 @@ def fit( elif self.dvals is not None: dvals = self.dvals else: - dvals = default_dose_grid(all_treated_doses) + dvals = default_dose_grid(modelled_doses) - # Build B-spline knots from all treated doses (unused on the discrete - # branch, but harmless to construct). + # Build B-spline knots from the modelled treated doses (excludes the d_L + # reference group under lowest_dose; unused on the discrete branch, but + # harmless to construct). knots, degree = build_bspline_basis( - all_treated_doses, degree=self.degree, num_knots=self.num_knots + modelled_doses, degree=self.degree, num_knots=self.num_knots ) # 3. Iterate over (g,t) cells @@ -976,6 +1103,20 @@ def fit( for gt, r in gt_results.items(): clean_gt[gt] = {k: v for k, v in r.items() if not k.startswith("_")} + # Unit-count metadata. Under lowest_dose the d_L group is the control / + # reference (not treated), so report the modelled-treated count (dose>d_L) + # and the reference-group size; reference_dose carries d_L. On the other + # paths the counts and reference_dose are unchanged (byte-stable). + if lowest_dose is not None: + _ud = treated_unit_doses.to_numpy(dtype=float) + n_treated_units_out = int(np.sum(_ud - lowest_dose > SATURATED_TOL)) + n_control_units_out = int(np.sum(np.abs(_ud - lowest_dose) <= SATURATED_TOL)) + reference_dose_out: Optional[float] = lowest_dose + else: + n_treated_units_out = int((unit_cohort > 0).sum()) + n_control_units_out = n_control + reference_dose_out = None + return ContinuousDiDResults( dose_response_att=dose_response_att, dose_response_acrt=dose_response_acrt, @@ -994,8 +1135,9 @@ def fit( groups=treatment_groups, time_periods=time_periods, n_obs=len(df), - n_treated_units=int((unit_cohort > 0).sum()), - n_control_units=n_control, + n_treated_units=n_treated_units_out, + n_control_units=n_control_units_out, + reference_dose=reference_dose_out, alpha=self.alpha, control_group=self.control_group, covariates=self.covariates, @@ -1347,6 +1489,7 @@ def _compute_dose_response_gt( dose_vector = precomp["dose_vector"] never_treated_mask = precomp["never_treated_mask"] time_periods = precomp["time_periods"] + lowest_dose = precomp.get("lowest_dose") # d_L reference (Remark 3.1) or None # Base period selection is_post = t >= g - self.anticipation @@ -1369,20 +1512,35 @@ def _compute_dose_response_gt( col_t = period_to_col[t] col_base = period_to_col[base_t] - # Treated units: first_treat == g and dose > 0 - treated_mask = (unit_cohorts == g) & (dose_vector > 0) + # Treated units: first_treat == g and dose > 0. Under lowest_dose + # (Remark 3.1) the lowest-dose group d_L is the comparison, so treated = + # doses strictly above d_L and the d_L group is the control. + if lowest_dose is not None: + treated_mask = (unit_cohorts == g) & (dose_vector - lowest_dose > SATURATED_TOL) + else: + treated_mask = (unit_cohorts == g) & (dose_vector > 0) n_treated = int(np.sum(treated_mask)) if n_treated == 0: return None - # Control units + # Control units (fail-closed dispatch so a new control_group value can + # never silently fall through to a wrong comparison group). if self.control_group == "never_treated": control_mask = never_treated_mask - else: - # Not-yet-treated: never-treated + first_treat > t + elif self.control_group == "not_yet_treated": control_mask = never_treated_mask | ( (unit_cohorts > t + self.anticipation) & (unit_cohorts != g) ) + elif self.control_group == "lowest_dose": + # Within-cohort lowest-dose group (single-cohort is enforced upstream, + # so this equals the pooled d_L group). d_L units are themselves + # treated at dose d_L, so subtracting their ΔY removes ATT(d_L) -> + # ATT(d)-ATT(d_L). + control_mask = (unit_cohorts == g) & ( + np.abs(dose_vector - lowest_dose) <= SATURATED_TOL + ) + else: # pragma: no cover - guarded by _validate_constrained_params + raise ValueError(f"Unhandled control_group: {self.control_group!r}") n_control = int(np.sum(control_mask)) if n_control == 0: warnings.warn( @@ -1459,8 +1617,13 @@ def _compute_dose_response_gt( def _design(z: np.ndarray) -> np.ndarray: return saturated_design_matrix(z, levels) + # Under lowest_dose the omitted reference is d_L (ATT(d_L)=0), so the + # backward-difference ACRT at the lowest modelled level references d_L + # (ACRT(d_1)=ATT(d_1)/(d_1-d_L)); base=0.0 otherwise (backward-to-zero). + _deriv_base = lowest_dose if lowest_dose is not None else 0.0 + def _deriv(z: np.ndarray) -> np.ndarray: - return saturated_derivative_design_matrix(z, levels) + return saturated_derivative_design_matrix(z, levels, base=_deriv_base) else: diff --git a/diff_diff/continuous_did_bspline.py b/diff_diff/continuous_did_bspline.py index c80ba2d8..7e8c06a3 100644 --- a/diff_diff/continuous_did_bspline.py +++ b/diff_diff/continuous_did_bspline.py @@ -321,32 +321,43 @@ def saturated_design_matrix(x, levels, tol=SATURATED_TOL): return B -def saturated_derivative_design_matrix(x, levels, tol=SATURATED_TOL): +def saturated_derivative_design_matrix(x, levels, tol=SATURATED_TOL, base=0.0): """ Finite-difference derivative rows for the saturated (discrete) basis. ACRT for a discrete dose is the paper's backward difference of the level - effects (CGBS 2024 §3.2 / §4.1) on the grid ``{d_0 = 0, d_1, ..., d_J}``, - where ``d_0 = 0`` is the omitted (untreated) category with ``ATT(0) = 0``: + effects (CGBS 2024 §3.2 / §4.1) on the grid ``{base, d_1, ..., d_J}``, where + ``base`` is the omitted reference category with ``ATT(base) = 0``: ``ACRT(d_j) = [ATT(d_j) - ATT(d_{j-1})] / (d_j - d_{j-1})``. At the lowest - positive level this references the zero-dose baseline, - ``ACRT(d_1) = [ATT(d_1) - 0] / (d_1 - 0) = ATT(d_1) / d_1`` — so a single - positive dose (``J = 1``, e.g. binary ``D in {0, 1}``) gives - ``ACRT(d_1) = ATT(d_1) / d_1`` and, for ``d_1 = 1``, the documented binary - identity ``ACRT = ATT``. This is a linear operator ``L`` on ``beta``, and - ``acrt = L @ beta``. Only the lowest level's row references ``d_0 = 0`` (so - that row does NOT sum to 0); the ``j >= 2`` rows are ordinary adjacent - backward differences (rows sum to 0). Returns the ``L`` row for each ``x_i`` - at its dose level. Analogous to :func:`bspline_derivative_design_matrix`. + modelled level this references ``base``, + ``ACRT(d_1) = [ATT(d_1) - 0] / (d_1 - base) = ATT(d_1) / (d_1 - base)``. + + ``base = 0.0`` (default) is the untreated baseline ``d_0 = 0`` for the D=0 + control path: ``ACRT(d_1) = ATT(d_1) / d_1`` — so a single positive dose + (``J = 1``, e.g. binary ``D in {0, 1}``) gives ``ACRT(d_1) = ATT(d_1) / d_1`` + and, for ``d_1 = 1``, the documented binary identity ``ACRT = ATT``. + ``base = d_L`` is the lowest-dose reference for the Remark 3.1 + (``control_group="lowest_dose"``) path, where the omitted category is the + lowest-dose group and ``ATT(d_L) = 0``: ``ACRT(d_1) = ATT(d_1) / (d_1 - d_L)``. + + This is a linear operator ``L`` on ``beta``, and ``acrt = L @ beta``. Only + the lowest level's row references ``base`` (so that row does NOT sum to 0); + the ``j >= 2`` rows are ordinary adjacent backward differences (rows sum to + 0). Returns the ``L`` row for each ``x_i`` at its dose level. Analogous to + :func:`bspline_derivative_design_matrix`. Parameters ---------- x : array-like Evaluation points, shape ``(n,)``. levels : array-like - Distinct dose levels, shape ``(J,)``. + Distinct modelled dose levels (excluding ``base``), shape ``(J,)``. tol : float, default=:data:`SATURATED_TOL` Matching tolerance. + base : float, default=0.0 + The omitted reference dose (``ATT(base) = 0``). ``0.0`` for the D=0 + control path; the lowest dose ``d_L`` for the lowest-dose-as-control + (Remark 3.1) path. Returns ------- @@ -357,12 +368,12 @@ def saturated_derivative_design_matrix(x, levels, tol=SATURATED_TOL): idx = _match_levels(x, levels, tol) J = len(levels) L = np.zeros((J, J)) - # Row 0 (lowest positive dose d_1): backward difference to the zero-dose - # baseline d_0 = 0, ATT(0) = 0 -> ACRT(d_1) = ATT(d_1) / d_1. Rows j >= 1: - # ordinary adjacent backward differences between positive doses. + # Row 0 (lowest modelled dose d_1): backward difference to the reference + # ``base`` with ATT(base) = 0 -> ACRT(d_1) = ATT(d_1) / (d_1 - base). Rows + # j >= 1: ordinary adjacent backward differences between modelled doses. for j in range(J): if j == 0: - L[0, 0] = 1.0 / levels[0] + L[0, 0] = 1.0 / (levels[0] - base) else: h = levels[j] - levels[j - 1] L[j, j - 1] = -1.0 / h diff --git a/diff_diff/continuous_did_results.py b/diff_diff/continuous_did_results.py index 779d5ef2..e4a93446 100644 --- a/diff_diff/continuous_did_results.py +++ b/diff_diff/continuous_did_results.py @@ -151,6 +151,10 @@ class ContinuousDiDResults: # per-dose-level regression); the ``dose_grid`` holds the distinct dose # levels when discrete. treatment_type: str = "continuous" + # Lowest-dose reference d_L for ``control_group="lowest_dose"`` (Remark 3.1); + # the estimand is ``ATT(d) - ATT(d_L)`` and ``ATT(d_L) = 0`` by construction. + # ``None`` for the never/not-yet-treated (D=0 control) paths. + reference_dose: Optional[float] = None event_study_effects: Optional[Dict[int, Dict[str, Any]]] = field(default=None) # Survey design metadata (SurveyMetadata instance from diff_diff.survey) survey_metadata: Optional[Any] = field(default=None) @@ -234,6 +238,9 @@ def summary(self, alpha: Optional[float] = None) -> str: f"{'Control group:':<30} {self.control_group:>10}", f"{'Treatment type:':<30} {self.treatment_type:>10}", ] + # Lowest-dose reference (Remark 3.1): show d_L when it is the control. + if self.reference_dose is not None: + lines.append(f"{'Reference dose (d_L):':<30} {self.reference_dose:>10.4g}") # Basis metadata: B-spline degree/knots (continuous) or the number of # saturated dose levels (discrete). if self.treatment_type == "discrete": diff --git a/diff_diff/guides/llms-autonomous.txt b/diff_diff/guides/llms-autonomous.txt index 54fb7d2e..a5584a2c 100644 --- a/diff_diff/guides/llms-autonomous.txt +++ b/diff_diff/guides/llms-autonomous.txt @@ -124,8 +124,8 @@ view. Every field below appears as a top-level key in that dict. Required by `SyntheticDiD`, `SunAbraham`, `EfficientDiD` under both `assumption="PT-All"` and `assumption="PT-Post"` (unless `control_group="last_cohort"` is passed), and `ContinuousDiD` - (which requires `P(D=0) > 0` - Remark 3.1 lowest-dose-as-control - is not yet implemented). Preferred-but-optional by + under its default `control_group` (`P(D=0) > 0`; when `P(D=0) = 0` + pass `control_group="lowest_dose"` for Remark 3.1). Preferred-but-optional by `CallawaySantAnna` and `ChaisemartinDHaultfoeuille`. Always `False` for `"categorical"`. - **`has_always_treated: bool`** - at least one binary-treatment @@ -231,10 +231,10 @@ view. Every field below appears as a top-level key in that dict. this setup, so it cannot encode timing). Under the canonical setup, several facts on the dose column predict `ContinuousDiD.fit()` outcomes: - `has_never_treated == True` (proxy for `P(D=0) > 0` under both - `control_group="never_treated"` and - `control_group="not_yet_treated"`; Remark 3.1 - lowest-dose-as-control is not yet implemented); + `has_never_treated == True` (proxy for `P(D=0) > 0` under the + default `control_group="never_treated"` / + `control_group="not_yet_treated"`; when `False`, pass + `control_group="lowest_dose"` for Remark 3.1); `treatment_varies_within_unit == False` (the actual fit-time gate, matching `ContinuousDiD.fit()`'s `df.groupby(unit)[dose].nunique() > 1` rejection at line @@ -248,9 +248,10 @@ view. Every field below appears as a top-level key in that dict. `dose_min` over non-zero values is the smallest treated dose). When `has_never_treated == False` (no zero-dose controls but - all observed doses non-negative), `ContinuousDiD` as currently - implemented does not apply (Remark 3.1 lowest-dose-as-control - is not implemented). Routing alternatives that do not require + all observed doses non-negative), pass + `control_group="lowest_dose"` (Remark 3.1) if there is a mass + point at the minimum dose — the lowest-dose group becomes the + comparison. Otherwise routing alternatives that do not require `P(D=0) > 0`: `HeterogeneousAdoptionDiD` for graded-adoption designs (HAD's own contract requires non-negative dose, which this branch satisfies), or linear DiD with the treatment as a @@ -515,9 +516,12 @@ When `has_never_treated == False`: (post-treatment periods at/after that cohort's adoption are trimmed). Distinct from CallawaySantAnna's `not_yet_treated` option. -- `ContinuousDiD` requires zero-dose control units (`P(D=0) > 0`). - Remark 3.1 of the paper (lowest-dose-as-control) is not yet - implemented; `fit()` raises `ValueError` when no `D=0` units exist. +- `ContinuousDiD` requires zero-dose control units (`P(D=0) > 0`) + under the default `control_group`; `fit()` raises `ValueError` when + no `D=0` units exist. When `P(D=0) = 0`, pass + `control_group="lowest_dose"` (Remark 3.1): the lowest-dose group + becomes the comparison (needs a mass point at the minimum dose), + estimand `ATT(d) - ATT(d_L)`. - `CallawaySantAnna` - use `control_group="not_yet_treated"` to use not-yet-treated units as the control pool. - `ChaisemartinDHaultfoeuille` - constructs switchers vs. non-switchers @@ -574,10 +578,10 @@ When `treatment_type == "continuous"`: `fit()` outcomes (full discussion in the paragraph immediately below): (a) zero-dose control units must exist (`PanelProfile.has_never_treated == True`, proxying - `ContinuousDiD`'s `P(D=0) > 0` requirement under both - `control_group` options because Remark 3.1 lowest-dose-as-control - is not yet implemented); (b) dose must be time-invariant per unit - (rule out panels where + `ContinuousDiD`'s `P(D=0) > 0` requirement under the default + `control_group` options; when it is `False`, pass + `control_group="lowest_dose"` for Remark 3.1); (b) dose must be + time-invariant per unit (rule out panels where `PanelProfile.treatment_varies_within_unit == True`); (c) the panel must be balanced (`PanelProfile.is_balanced == True`); (d) no `duplicate_unit_time_rows` alert (the precompute path @@ -867,8 +871,8 @@ Reasoning chain: column). On the dose column alone, profile_panel exposes five facts that predict `fit()` outcomes under that canonical setup: `has_never_treated == True` (proxy for `P(D=0) > 0` - under both `control_group` options, since Remark 3.1 - lowest-dose-as-control is not yet implemented), + under the default `control_group` options; when `False`, pass + `control_group="lowest_dose"` for Remark 3.1), `treatment_varies_within_unit == False` (the actual fit-time gate matching `ContinuousDiD.fit()`'s `df.groupby(unit)[dose].nunique() > 1` rejection at line @@ -904,12 +908,12 @@ Reasoning chain: 4. Counter-example: had `has_never_treated == False` (every unit eventually treated, even if some pre-treatment rows have zero dose so `treatment_dose.has_zero_dose == True`), - `ContinuousDiD.fit()` would reject the panel under both - `control_group="never_treated"` and - `control_group="not_yet_treated"` because Remark 3.1 - lowest-dose-as-control is not yet implemented. On this branch - (no never-treated controls but doses still non-negative), - `HeterogeneousAdoptionDiD` IS a routing alternative for + `ContinuousDiD.fit()` would reject the panel under + `control_group="never_treated"` / `control_group="not_yet_treated"` + (both require `P(D=0) > 0`); pass `control_group="lowest_dose"` + (Remark 3.1) to use the lowest-dose group as the comparison. On + this branch (no never-treated controls but doses still + non-negative), `HeterogeneousAdoptionDiD` IS a routing alternative for graded-adoption designs, and linear DiD with the treatment as a continuous covariate is another; see §2 for the full routing discussion. diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index c588e1a8..9f921ec7 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -688,7 +688,7 @@ ContinuousDiD( degree: int = 3, # B-spline degree (3 = cubic) num_knots: int = 0, # Interior knots dvals: np.ndarray | None = None, # Custom dose evaluation grid - control_group: str = "never_treated", # "never_treated" or "not_yet_treated" + control_group: str = "never_treated", # "never_treated", "not_yet_treated", or "lowest_dose" (Remark 3.1, P(D=0)=0) anticipation: int = 0, base_period: str = "varying", # "varying" or "universal" alpha: float = 0.05, @@ -2367,7 +2367,7 @@ DIFF_DIFF_BACKEND=rust pytest # Force Rust (fail if unavailable) | Staggered treatment timing | `CallawaySantAnna`, `ImputationDiD`, or `SunAbraham` | | Few treated units / synthetic control | `SyntheticDiD` | | Interactive fixed effects / factor confounding | `TROP` | -| Continuous treatment intensity, per-dose ATT(d) / ACRT(d) (requires never-treated controls) | `ContinuousDiD` | +| Continuous treatment intensity, per-dose ATT(d) / ACRT(d) (needs zero-dose controls, or `control_group="lowest_dose"` for Remark 3.1 when P(D=0)=0) | `ContinuousDiD` | | Continuous treatment intensity, WAS at dose support boundary (compatible with universal rollout or small never-treated share) | `HeterogeneousAdoptionDiD` | | Two-criterion treatment, simultaneous (2x2x2 DDD) | `TripleDifference` | | Two-criterion treatment, staggered timing + eligibility | `StaggeredTripleDifference` | diff --git a/diff_diff/practitioner.py b/diff_diff/practitioner.py index 3a636e3c..54780379 100644 --- a/diff_diff/practitioner.py +++ b/diff_diff/practitioner.py @@ -1059,8 +1059,10 @@ def _handle_had(results: Any): "throughout, encoded as either first_treat == 0 or " "first_treat == inf; ContinuousDiD normalizes inf -> 0 " "internally), ContinuousDiD is the alternative — " - "different estimand, and ContinuousDiD's identification " - "requires never-treated controls. HAD itself remains " + "different estimand; ContinuousDiD's default " + "identification uses never-treated controls (or " + "control_group='lowest_dose' for Remark 3.1 when " + "P(D=0)=0). HAD itself remains " "valid even with a small share of never-treated units " "(paper compatibility; see REGISTRY § " "HeterogeneousAdoptionDiD edge cases — Garrett et al. " @@ -1238,7 +1240,9 @@ def _handle_had_event_study(results: Any): "either first_treat == 0 or first_treat == inf; " "ContinuousDiD normalizes inf -> 0 internally), " "ContinuousDiD is the alternative — different estimand, " - "requires never-treated. Two ContinuousDiD aggregation " + "uses never-treated by default (or " + "control_group='lowest_dose' for Remark 3.1 when " + "P(D=0)=0). Two ContinuousDiD aggregation " "surfaces are relevant and distinct: `aggregate='dose'` " "(the default) produces the per-dose ATT(d) / ACRT(d) " "curves on `results.dose_response_att` / " diff --git a/diff_diff/profile.py b/diff_diff/profile.py index 62c6b4b8..5b0548fa 100644 --- a/diff_diff/profile.py +++ b/diff_diff/profile.py @@ -81,13 +81,14 @@ class TreatmentDoseShape: 1. ``PanelProfile.has_never_treated == True`` (some unit has dose 0 in every period). Predicts the estimator's - ``P(D=0) > 0`` requirement under both - ``control_group="never_treated"`` and - ``control_group="not_yet_treated"`` (Remark 3.1 - lowest-dose-as-control not yet implemented), because the - canonical setup ties ``first_treat == 0`` to ``D_i == 0``. - Failure means no never-treated controls exist on the dose - column; see routing notes below. + ``P(D=0) > 0`` requirement under the default + ``control_group="never_treated"`` / + ``control_group="not_yet_treated"`` (the canonical setup ties + ``first_treat == 0`` to ``D_i == 0``). When it is ``False``, + ``control_group="lowest_dose"`` (Remark 3.1) is the route — + the lowest-dose group becomes the comparison (needs a mass + point at the minimum dose) — so failure of (1) no longer rules + ``ContinuousDiD`` out; see routing notes below. 2. ``PanelProfile.treatment_varies_within_unit == False`` (per-unit full-path dose constancy on the dose column). This IS the actual fit-time gate, matching @@ -114,12 +115,14 @@ class TreatmentDoseShape: Routing alternatives when (1) or (5) fails: - When (1) fails (no never-treated controls but all observed - doses non-negative): ``ContinuousDiD`` does not apply (Remark - 3.1 lowest-dose-as-control is not implemented). - ``HeterogeneousAdoptionDiD`` IS a candidate for graded-adoption - designs (HAD's contract requires non-negative dose, satisfied - here); linear DiD with the treatment as a continuous covariate - is another. + doses non-negative): use ``control_group="lowest_dose"`` + (Remark 3.1) if there is a mass point at the minimum dose + (``>= 2`` units at ``d_L``) — the lowest-dose group becomes the + comparison and the estimand is ``ATT(d) - ATT(d_L)``. + Otherwise ``HeterogeneousAdoptionDiD`` IS a candidate for + graded-adoption designs (HAD's contract requires non-negative + dose, satisfied here); linear DiD with the treatment as a + continuous covariate is another. - When (5) fails (negative treated doses): ``HeterogeneousAdoptionDiD`` is **not** a fallback either — HAD raises on negative post-period dose (``had.py:1450-1459``, diff --git a/docs/api/continuous_did.rst b/docs/api/continuous_did.rst index fc3b758b..ad442ff0 100644 --- a/docs/api/continuous_did.rst +++ b/docs/api/continuous_did.rst @@ -158,5 +158,5 @@ Comparison with CallawaySantAnna - Multiplier bootstrap (optional) - Multiplier bootstrap (optional) * - Control group - - never_treated / not_yet_treated + - never_treated / not_yet_treated / lowest_dose - never_treated / not_yet_treated diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 39a80873..8b66c66c 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -961,6 +961,7 @@ See `docs/methodology/continuous-did.md` Section 4 for full details. - `ATT^{loc} = E[ATT(D|D) | D > 0] = E[Delta Y | D > 0] - E[Delta Y | D = 0]` — binarized ATT (PT); equals `ATT^{glob}` under SPT - `ATT^{glob} = E[ATT(D) | D > 0]` — global average dose-response level (SPT required) - `ACRT^{glob} = E[ACRT(D_i) | D > 0]` — plug-in average marginal effect (SPT required) +- **Lowest-dose-as-control (Remark 3.1, `control_group="lowest_dose"`):** when `P(D=0) = 0` (no untreated group) the lowest-dose group `d_L` is the comparison and the targets are rebased to it — `ATT(d) - ATT(d_L)` (with `ATT(d_L) = 0` by construction, the omitted reference), `ACRT(d)` backward-differenced to `d_L` (`ACRT(d_1) = ATT(d_1)/(d_1 - d_L)`), and `ATT^{glob} = E[ΔY | D > d_L] - E[ΔY | D = d_L]`. Requires a genuine lowest-dose group (`P(D = d_L) > 0`). See Note #7. **Estimation via B-spline OLS:** 1. Compute `Delta_tilde_Y = (Y_t - Y_{t-1})_treated - mean((Y_t - Y_{t-1})_control)` @@ -1002,7 +1003,7 @@ there (analytical `ACRT(d_1)` SE matches the bootstrap). Otherwise the ### Edge Cases -- **No untreated group**: Remark 3.1 (lowest-dose-as-control) not implemented; requires P(D=0) > 0. +- **No untreated group**: `control_group="lowest_dose"` implements Remark 3.1 — the lowest-dose group `d_L` becomes the comparison (estimand `ATT(d) - ATT(d_L)`); requires a genuine lowest-dose group (`P(D=d_L) > 0`, `>= 2` units at `d_L`) and no never-treated units. Single-cohort only (multi-cohort and `covariates=` raise; see Note #7). The default `never_treated`/`not_yet_treated` paths still require `P(D=0) > 0`. - **Discrete treatment**: `treatment_type="discrete"` fits the saturated per-dose-level regression (CGBS 2024 Eq. 4.1; see Note #6). On the default `treatment_type="continuous"` path an integer-valued dose is detected and warns, pointing to `treatment_type="discrete"`. - **All-same dose**: B-spline basis collapses; ACRT(d) = 0 everywhere. - **Rank deficiency**: When n_treated <= n_basis, cell is skipped. @@ -1034,6 +1035,7 @@ labels.* 4. **Note:** Zero-`first_treat` rows with nonzero `dose` are force-zeroed with `UserWarning` reporting the affected row count (axis-E silent-coercion). No R correspondence; `contdid` v0.1.0 has the same `first_treat = 0` → `D = 0` invariant but silently coerces without a warning. Cross-references the § Implementation Checklist `**Note:**` below and `METHODOLOGY_REVIEW.md` § ContinuousDiD Deviations #4. 5. **Note (covariate support — library extension beyond `contdid` v0.1.0):** `covariates=` with `estimation_method ∈ {reg, dr}` adds conditional-parallel-trends adjustment. This is a **library extension**: `contdid` v0.1.0 hard-stops on any covariate (`stop("covariates not currently supported…")`), so there is **no external R anchor for the covariate-adjusted dose *curve***. Validation instead: (a) the **scalar `overall_att` + SE** map *exactly* onto `DRDID::reg_did_panel` (reg) / `DRDID::drdid_panel` (dr) — a tight (~1e-8) component anchor, skip-guarded since DRDID is not in CI (`tests/test_methodology_continuous_did.py::TestCovariateReg`); (b) an **R-free NumPy reconstruction** of the reg/dr `att`+SE runs *in CI* at p≥2 (`test_dr_reg_numpy_crosscheck_p2`) — the guard the p=1 reduction cannot provide (at p=1 the intercept-only propensity is constant, so `eta_cont ≡ 0` and dr collapses to reg); (c) DGP recovery + MC coverage (reg 96%, dr 95%). **`ipw` restricted:** `estimation_method="ipw"` with covariates raises `NotImplementedError` — pure IPW's covariate adjustment is a single scalar (a propensity-reweighted control mean) that shifts only the ATT(d) level and leaves `ACRT(d)` identical to the unconditional fit, so it cannot adjust the dose-response *shape*. **Deviations from DRDID:** unit weights are 1 (unweighted; `covariates=` + `survey_design=` raises `NotImplementedError`, deferred); propensity trimming uses clip semantics (`pscore_trim`) rather than DRDID's drop-trimming — identical on moderate-overlap data (the anchor regime), diverging only at extreme propensities. **Fail-closed policies (no-silent-failures):** (i) missing/non-finite covariate values raise `ValueError` up front — a per-cell fallback to unconditional estimation would silently mix conditional-PT and unconditional-PT cells in the aggregate; (ii) `dr` propensity-estimation failure raises by default (`pscore_fallback="error"`) so a `dr` fit never silently degrades to a non-DR estimate — `pscore_fallback="unconditional"` opts into the graceful (warned, reg-like) fallback. **`treatment_type="discrete"` composability:** on the intercept-free saturated basis the constant augmentation `eta_cont` shifts every `beta_j` equally (indicator partition of unity), so it cancels in the adjacent differences and reg/dr share `ACRT(d_j)` point AND SE for `j >= 2` — but the lowest-dose ACRT references the fixed baseline `ATT(0) = 0` (backward-to-zero), so reg/dr **differ at `ACRT(d_1)` by `eta_cont/d_1`** (the dr IF carries the augmentation variance there). See Note #6. Cross-references `docs/methodology/continuous-did.md` § Covariates. 6. **Note (discrete-treatment saturated regression — library extension beyond `contdid` v0.1.0):** `treatment_type="discrete"` estimates the dose-response by a **saturated regression** (CGBS 2024 Eq. 4.1) — one indicator per distinct dose level, so `beta_j = mean_{D=d_j}(ΔY − control) = ATT(d_j)` (a per-level 2×2 DiD) — instead of the B-spline sieve. `ACRT(d_j)` is the paper's **backward difference** on the grid `{d_0 = 0, d_1, …, d_J}` (Eq. 4.1 makes `d_0 = 0` the omitted category with `ATT(0) = 0`): `ACRT(d_j) = [ATT(d_j) − ATT(d_{j-1})]/(d_j − d_{j-1})` for `j ≥ 2`, and at the lowest positive level it references the zero-dose baseline, `ACRT(d_1) = [ATT(d_1) − 0]/(d_1 − 0) = ATT(d_1)/d_1`. So a single positive dose (`J = 1`, e.g. binary `D ∈ {0,1}`) yields `ACRT(d_1) = ATT(d_1)/d_1`, and for `d_1 = 1` the documented binary identity `ACRT = ATT` holds exactly. This is a **library extension**: `contdid` v0.1.0 accepts `treatment_type` in its signature but **does not implement the discrete path** (documented "Discrete treatment not yet implemented"), so there is **no external R anchor**. It is instead an *exact* basis swap of the B-spline design/evaluation/derivative trio for an indicator/identity/finite-difference trio; every downstream quantity is linear in `beta`, so the analytical-SE / multiplier-bootstrap / covariate (reg,dr) / survey machinery is reused unchanged and reduces *analytically* to the per-level 2×2 DiD (`bread @ psi_bar = ones(J)`; the common control mean cancels in the `j ≥ 2` adjacent differences whose `L`-rows sum to 0). **reg vs dr:** the constant DR augmentation `η̄_cont` cancels in the `j ≥ 2` differences, so `ACRT(d_j)` point AND SE are identical for `reg`/`dr` there; but `ACRT(d_1) = ATT(d_1)/d_1` references the fixed baseline `ATT(0) = 0` (not shifted by `η̄_cont`), so `reg` and `dr` genuinely **differ at `ACRT(d_1)` by `η̄_cont/d_1`** (and correspondingly in `ACRT^glob` via the `d_1` mass) — the dr influence function carries the augmentation variance at `d_1` (validated: analytical `ACRT(d_1)` SE matches the multiplier bootstrap). Validation (R-free, in CI): exact hand-calc of `ATT(d_j)`/`ACRT`/`overall_att` and the analytical SE against a direct per-level 2×2 reconstruction (`~1e-12`/`~1e-10`), DGP recovery, and MC coverage for analytical + bootstrap (`tests/test_methodology_continuous_did.py::TestDiscreteSaturated`, `tests/test_continuous_did.py::TestDiscreteSaturatedAPI`). **Fail-closed policies (no-silent-failures):** (i) multi-cohort fits with **heterogeneous dose support** across cohorts raise `NotImplementedError` — an absent global level yields a dropped zero column (`att_d[level]=0`) that the plain-sum dose aggregation would bias toward zero (support-aware aggregation is deferred; single-cohort, 2-period, and shared-support multi-cohort are supported); (ii) a requested `dvals` value that is not an observed dose level raises `ValueError` (a saturated model cannot be evaluated off-support); (iii) an over-parameterized fit (`< 2` treated units per level, or `J > n_treated/2`) warns (degenerate per-level SE); (iv) with `survey_design=`, any dose level with **zero effective treated mass in a `(g,t)` cell** raises `ValueError` — a per-cell check (not just the global positive-weight check), so a level that survey/subpopulation weights zero out for one cohort while another cohort keeps it cannot silently drop to a zero-coefficient saturated column. Cross-references `docs/methodology/continuous-did.md` § 5.1. +7. **Note (lowest-dose-as-control, Remark 3.1 — library extension beyond `contdid` v0.1.0):** `control_group="lowest_dose"` implements CGBS 2024 Remark 3.1 for settings with no untreated group (`P(D=0) = 0`): the lowest-dose group `d_L` becomes the comparison and the estimand is `ATT(d) − ATT(d_L)` (SPT), with `ATT(d_L) = 0` the omitted reference. Mechanically it is a **control-group swap** — the D=0 control pool is replaced by the `d_L` group; the entire linear influence-function / bootstrap / event-study / survey machinery is control-group-generic and reused unchanged (`ee_control` already carries the reference-group variance, so **no new SE plumbing**). On the discrete saturated basis the backward-difference operator's reference shifts from `0` to `d_L` (`ACRT(d_1) = ATT(d_1)/(d_1 − d_L)`); on the continuous B-spline path the reference shifts only `μ_0` (the level), leaving `ACRT = spline'` unchanged. `contdid` v0.1.0 does **not** implement Remark 3.1, so there is **no external R anchor**; validation (R-free, in CI): an **exact `d_L → 0` equivalence** anchor (relabelling a `never_treated` panel's D=0 group as a tiny common dose `d_L = ε` reproduces the `never_treated` ATT and SE exactly, for any ε), a discrete hand-calc of `ATT(d)−ATT(d_L)`/`ACRT`/`overall_att`/`overall_acrt` and the per-level 2×2 SE (`~1e-10`), continuous mass-point DGP recovery, analytical-vs-bootstrap SE agreement, a pre-period placebo, and MC coverage (`tests/test_methodology_continuous_did.py::TestLowestDose`, `tests/test_continuous_did.py::TestLowestDoseAPI`). The continuous path requires a genuine **mass point** at the minimum dose (`>= 2` units at `d_L`, i.e. `P(D=d_L) > 0`) — the Remark 3.1 identification condition; a singleton minimum fails closed. **Fail-closed policies (no-silent-failures):** (i) never-treated units present with `lowest_dose` → `ValueError` (they would be silently dropped); (ii) singleton `d_L` (no mass point) → `ValueError`; (iii) no treated dose above `d_L` → `ValueError`; (iv) user `dvals ≤ d_L` → `ValueError` (`d_L` is the omitted reference); (v) survey/subpopulation weighting that leaves the `d_L` group with `< 2` positive-weight units → `ValueError` (a single positive-weight reference unit gives `ee_control = 0`, i.e. zero control-side variance — the effective-`>= 2` analogue of the raw mass-point guard, applied after weighting); (vi) a boundary gap `d_1 − d_L` that is a tiny fraction of the dose range warns (huge boundary ACRT/SE). **Deferred (fail-closed `NotImplementedError` + TODO):** multi-cohort `lowest_dose` (needs a within-cohort reference + support-aware cross-cohort aggregation) and `covariates=` × `lowest_dose` (conditional-PT-relative-to-`d_L` estimand). Cross-references `docs/methodology/continuous-did.md` § 5.6. ### Implementation Checklist @@ -1043,9 +1045,9 @@ labels.* - [x] Multiplier bootstrap for inference - [x] Analytical SEs via influence functions - [x] Equation verification tests (linear, quadratic, multi-period) -- [x] Covariate support (reg / dr) — conditional parallel trends. **ipw restricted** (see Note below); survey × covariate deferred; Remark 3.1 still deferred. +- [x] Covariate support (reg / dr) — conditional parallel trends. **ipw restricted** (see Note below); survey × covariate deferred; covariates × lowest_dose deferred. - [x] Discrete treatment saturated regression (`treatment_type="discrete"`) — CGBS 2024 Eq. 4.1; per-level ATT(d_j) + backward-difference ACRT on `{0, d_1, …, d_J}` (`ACRT(d_1) = ATT(d_1)/d_1`; binary → ACRT = ATT). Multi-cohort heterogeneous dose support deferred (see Note #6). -- [ ] Lowest-dose-as-control (Remark 3.1) +- [x] Lowest-dose-as-control (Remark 3.1, `control_group="lowest_dose"`) — estimand `ATT(d) - ATT(d_L)` for `P(D=0)=0`; discrete + continuous (mass-point) paths; single-cohort. Multi-cohort and covariates × lowest_dose deferred (see Note #7). - [x] Survey design support (Phase 3): weighted B-spline OLS, TSL on influence functions; bootstrap+survey supported (Phase 6) - **Note:** ContinuousDiD bootstrap with survey weights supported (Phase 6) via PSU-level multiplier weights - **Note:** The R-style convention of coding never-treated units as `first_treat=inf` is still accepted and normalized to `first_treat=0` internally, but the estimator now emits a `UserWarning` reporting the row count so the silent recategorization is surfaced (axis-E silent coercion under the Phase 2 audit). Only `+inf` is recoded (matching the R convention). Any **negative** `first_treat` value (including `-inf`) raises `ValueError` with the row count, since such units would otherwise silently fall out of both the treated (`g > 0`) and never-treated (`g == 0`) masks. Pass `0` directly for never-treated units to avoid the warning. diff --git a/docs/methodology/continuous-did.md b/docs/methodology/continuous-did.md index 33d73ee7..aac1fefd 100644 --- a/docs/methodology/continuous-did.md +++ b/docs/methodology/continuous-did.md @@ -293,6 +293,16 @@ When P(D=0) = 0 (all units receive some treatment), use the lowest dose group d_ as comparison. Under PT, this recovers ATT(d|d) - ATT(d_L|d_L). Under SPT, recovers ATT(d) - ATT(d_L). +**Implemented** via `control_group="lowest_dose"` (both `treatment_type="discrete"` +and the continuous B-spline path). Mechanically a control-group swap: the D=0 +control pool is replaced by the d_L group, so `mu_0 = mean(ΔY | D = d_L)` and the +per-cell estimand becomes `ATT(d) - ATT(d_L)` (with `ATT(d_L) = 0` the omitted +reference). The continuous path requires a genuine mass point at d_L +(`P(D=d_L) > 0`, the Remark 3.1 identification condition); a singleton minimum +fails closed. Single-cohort only in v1 — multi-cohort (which needs a within-cohort +reference and support-aware cross-cohort aggregation) and `covariates=` × +`lowest_dose` (conditional PT relative to d_L) raise `NotImplementedError`. + --- ## 6. Inference @@ -472,3 +482,4 @@ can be meaningfully aggregated. - TWFE decomposition diagnostics *Discrete treatment (saturated regression) is now implemented — see § 5.1.* +*Lowest-dose-as-control (Remark 3.1) is now implemented via `control_group="lowest_dose"` — see § 5.6.* diff --git a/docs/practitioner_decision_tree.rst b/docs/practitioner_decision_tree.rst index 977080fc..5253ea1a 100644 --- a/docs/practitioner_decision_tree.rst +++ b/docs/practitioner_decision_tree.rst @@ -279,11 +279,17 @@ Universal Rollout (No Untreated Markets) ---------------------------------------- **Your situation:** Every market got the campaign at the same time - there is no -holdout group - but spending levels varied across markets. ``ContinuousDiD`` cannot -help here because it requires an untreated comparison group; standard DiD has no -control to anchor the contrast. - -**Recommended method:** :class:`~diff_diff.HeterogeneousAdoptionDiD` +holdout group - but spending levels varied across markets. ``ContinuousDiD``'s +default control groups require an untreated (zero-dose) comparison group; standard +DiD likewise has no control to anchor the contrast. The choice hinges on whether +the lowest spending level forms a genuine group (a mass point, ``>= 2`` markets at +the minimum). + +**Recommended method:** :class:`~diff_diff.ContinuousDiD` with +``control_group="lowest_dose"`` **when a lowest-dose group exists** and you want the +per-dose curve relative to that anchor (CGBS Remark 3.1, estimand +``ATT(d) - ATT(d_L)``); otherwise :class:`~diff_diff.HeterogeneousAdoptionDiD` for +the Weighted Average Slope target below. This estimator implements de Chaisemartin, Ciccia, D'Haultfoeuille and Knau (2026) and resolves to one of two estimands depending on whether the smallest-dose diff --git a/tests/test_continuous_did.py b/tests/test_continuous_did.py index 590b8661..bfe80e9b 100644 --- a/tests/test_continuous_did.py +++ b/tests/test_continuous_did.py @@ -2,6 +2,8 @@ Unit and integration tests for ContinuousDiD estimator. """ +import warnings + import numpy as np import pandas as pd import pytest @@ -1952,3 +1954,197 @@ def test_continuous_default_matches_explicit(self): r_default.dose_response_att.effects, r_explicit.dose_response_att.effects ) np.testing.assert_allclose(r_default.overall_att, r_explicit.overall_att) + + +class TestLowestDoseAPI: + """API, guards, and metadata for control_group='lowest_dose' (Remark 3.1).""" + + def _no_d0(self, effects=None, **kw): + """Discrete panel with NO never-treated units (P(D=0)=0).""" + return _discrete_panel(effects or {1.0: 0.5, 2.0: 1.5, 4.0: 2.5}, n_control=0, **kw) + + def test_lowest_dose_valid_in_params_transactional(self): + est = ContinuousDiD() + est.set_params(control_group="lowest_dose") + assert est.control_group == "lowest_dose" + assert est.get_params()["control_group"] == "lowest_dose" + # Transactional: an invalid update leaves config unmutated. + with pytest.raises(ValueError): + est.set_params(control_group="bogus") + assert est.control_group == "lowest_dose" + + def test_never_treated_present_raises(self): + """lowest_dose with never-treated units present fails closed (no silent drop).""" + df = _discrete_panel({1.0: 0.5, 2.0: 1.5}, n_control=40) # has D=0 units + est = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete") + with pytest.raises(ValueError, match="never-treated"): + est.fit(df, **_DKW) + + def test_singleton_dL_raises(self): + """A singleton minimum dose is not a lowest-dose group -> ValueError.""" + df = self._no_d0({2.0: 1.0, 4.0: 2.0}) + # Add a single unit at dose 1.0 (a singleton minimum). + extra = df[df["unit"] == df["unit"].iloc[0]].copy() + extra["unit"] = int(df["unit"].max()) + 1 + extra["dose"] = 1.0 + df = pd.concat([df, extra], ignore_index=True) + est = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete") + with pytest.raises(ValueError, match="lowest-dose"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + est.fit(df, **_DKW) + + def test_multi_cohort_raises(self): + """Multi-cohort lowest_dose is deferred -> NotImplementedError.""" + df = _discrete_panel({1.0: 0.5, 2.0: 1.5}, n_control=0, cohorts=(1, 2)) + est = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete") + with pytest.raises(NotImplementedError, match="multiple treatment cohorts"): + est.fit(df, **_DKW) + + def test_covariates_raises_at_init(self): + """covariates + lowest_dose is deferred -> NotImplementedError (config-level).""" + with pytest.raises(NotImplementedError, match="covariates"): + ContinuousDiD(control_group="lowest_dose", covariates=["x1"]) + + def test_dvals_below_dL_raises(self): + """User dvals at/below the reference d_L are rejected on both paths.""" + df = self._no_d0() + est = ContinuousDiD( + control_group="lowest_dose", treatment_type="discrete", dvals=[1.0, 2.0] + ) + with pytest.raises(ValueError, match="reference dose"): + est.fit(df, **_DKW) + + def test_continuous_no_mass_point_raises(self): + """Continuous dose with a singleton minimum (no mass point) -> ValueError.""" + rng = np.random.default_rng(3) + rows, uid = [], 0 + for d in np.linspace(1.0, 5.0, 80): # all distinct -> singleton min + base = rng.normal(0, 1) + for p in (0, 1): + rows.append((uid, p, base + 0.3 * p + rng.normal(0, 0.3), 1, float(d))) + uid += 1 + df = pd.DataFrame(rows, columns=["unit", "period", "outcome", "first_treat", "dose"]) + est = ContinuousDiD(control_group="lowest_dose", treatment_type="continuous", degree=2) + with pytest.raises(ValueError, match="mass point"): + est.fit( + df, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + dose="dose", + aggregate="dose", + ) + + def test_continuous_one_dose_above_warns(self): + """Exactly one modelled dose above d_L: warn + valid ATT level, ACRT=0 (M3).""" + rng = np.random.default_rng(4) + rows, uid = [], 0 + for d in [1.0] * 40 + [2.0] * 40: # mass point at 1, one dose above + base = rng.normal(0, 1) + for p in (0, 1): + y = base + 0.3 * p + (0.8 * (d - 1.0) if p >= 1 else 0.0) + rng.normal(0, 0.3) + rows.append((uid, p, y, 1, d)) + uid += 1 + df = pd.DataFrame(rows, columns=["unit", "period", "outcome", "first_treat", "dose"]) + est = ContinuousDiD(control_group="lowest_dose", treatment_type="continuous", degree=2) + with pytest.warns(UserWarning): + res = est.fit( + df, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + dose="dose", + aggregate="dose", + ) + assert np.all(np.isfinite(res.dose_response_att.se)) + assert np.allclose(res.dose_response_acrt.effects, 0.0) + + def test_threshold_boundary_no_phantom_level(self): + """A unit at d_L + 5e-10 clusters into the control group (no phantom level).""" + df = self._no_d0() # doses {1, 2, 4}, d_L = 1 + extra = df[df["unit"] == df["unit"].iloc[0]].copy() + extra["unit"] = int(df["unit"].max()) + 1 + extra["dose"] = 1.0 + 5e-10 # within SATURATED_TOL of d_L + df = pd.concat([df, extra], ignore_index=True) + res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit(df, **_DKW) + # No spurious level at ~1.0; modelled grid is exactly {2, 4}. + assert np.allclose(res.dose_response_att.dose_grid, [2.0, 4.0]) + + def test_metadata_and_summary(self): + df = self._no_d0() + res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit(df, **_DKW) + assert res.reference_dose == 1.0 + assert res.n_control_units == 40 # the d_L group (n_per=40) + assert res.n_treated_units == 80 # doses 2 and 4, 40 each + assert res.control_group == "lowest_dose" + s = res.summary() + assert "Reference dose" in s and "lowest_dose" in s + + def test_reference_dose_none_on_other_paths(self): + """reference_dose is None for never/not-yet-treated (byte-stable metadata).""" + df = _discrete_panel({1.0: 0.5, 2.0: 1.5}, n_control=40) + res = ContinuousDiD(control_group="never_treated", treatment_type="discrete").fit( + df, **_DKW + ) + assert res.reference_dose is None + + def test_reworded_not_yet_treated_error_points_at_lowest_dose(self): + """The no-D=0 not_yet_treated error still raises and points at lowest_dose (L5).""" + df = self._no_d0() + est = ContinuousDiD(control_group="not_yet_treated", treatment_type="discrete") + with pytest.raises(ValueError, match="lowest_dose"): + est.fit(df, **_DKW) + + def test_event_study_lowest_dose(self): + """Event-study aggregation composes with lowest_dose (the ES IF swaps control).""" + df = _discrete_panel( + {1.0: 0.5, 2.0: 1.5}, n_control=0, cohorts=(2,), n_periods=4, noise=0.0, seed=6 + ) + res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( + df, "outcome", "unit", "period", "first_treat", "dose", aggregate="eventstudy" + ) + assert res.event_study_effects is not None + # Pre-period event bins (e < 0) difference out to ~0 (both groups untreated). + pre = {e: v for e, v in res.event_study_effects.items() if e < 0} + assert pre and all(abs(v["effect"]) < 1e-9 for v in pre.values()) + + def test_survey_zeroed_dL_group_raises(self): + """A survey design zeroing the entire d_L reference group fails closed.""" + from diff_diff import SurveyDesign + + df = self._no_d0() # doses {1, 2, 4}, d_L = 1 + # Zero the survey weight of every d_L unit -> no reference group remains. + df["wt"] = np.where(np.abs(df["dose"] - 1.0) <= 1e-9, 0.0, 1.0) + est = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete") + with pytest.raises(ValueError, match="positive-weight unit"): + est.fit(df, survey_design=SurveyDesign(weights="wt"), **_DKW) + + def test_survey_effective_singleton_dL_raises(self): + """Survey weights leaving only ONE positive-weight d_L unit fail closed. + + The raw >= 2 guard runs before weighting; a subpopulation keeping a + single positive-weight d_L unit gives zero reference-side variance + (its ee_control = w*(dY - mu_0) = 0), so require >= 2 effective units. + """ + from diff_diff import SurveyDesign + + df = self._no_d0() # doses {1, 2, 4}, d_L = 1, >= 2 raw units at d_L + is_dL = np.abs(df["dose"] - 1.0) <= 1e-9 + keep = sorted(df.loc[is_dL, "unit"].unique())[0] # keep exactly one d_L unit + df["wt"] = np.where(is_dL & (df["unit"] != keep), 0.0, 1.0) + est = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete") + with pytest.raises(ValueError, match="positive-weight unit"): + est.fit(df, survey_design=SurveyDesign(weights="wt"), **_DKW) + + def test_idempotent_refit(self): + """fit() does not mutate config; clone + refit reproduces the fit.""" + df = self._no_d0() + est = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete") + cfg = est.get_params() + r1 = est.fit(df, **_DKW) + assert est.get_params() == cfg # config unchanged by fit + r2 = est.fit(df, **_DKW) + np.testing.assert_allclose(r1.dose_response_att.effects, r2.dose_response_att.effects) diff --git a/tests/test_methodology_continuous_did.py b/tests/test_methodology_continuous_did.py index 3b9e11a7..1b2e8ae2 100644 --- a/tests/test_methodology_continuous_did.py +++ b/tests/test_methodology_continuous_did.py @@ -9,6 +9,7 @@ import os import subprocess import tempfile +import warnings import numpy as np import pandas as pd @@ -1206,6 +1207,49 @@ def _hand_calc_discrete(df, levels): return att, acrt, overall, se +def _hand_calc_lowest_dose(df, all_levels): + """Independent NumPy lowest-dose-as-control (Remark 3.1) recomputation. + + ``all_levels`` are the sorted distinct positive doses; ``d_L = all_levels[0]`` + is the reference (control), modelled levels are ``all_levels[1:]``. Returns + ``(modelled_levels, att, acrt, overall, overall_acrt, se)`` where + ``att[j] = mean(dY | d_j) - mean(dY | d_L)`` (a per-level 2x2 DiD vs the + lowest-dose group), ACRT is the backward difference with base ``d_L``, and + ``se`` is the per-level 2x2 SE with the ``d_L`` group as control. + """ + periods = sorted(df["period"].unique()) + p0, p1 = periods[0], periods[-1] + wide = df.pivot(index="unit", columns="period", values="outcome") + dy = (wide[p1] - wide[p0]).to_numpy() + udose = df.groupby("unit")["dose"].first().to_numpy() + d_L = all_levels[0] + modelled = all_levels[1:] + cm = np.abs(udose - d_L) <= 1e-9 # d_L group = control + mu0 = dy[cm].mean() + att = np.array([dy[np.abs(udose - d) <= 1e-9].mean() - mu0 for d in modelled]) + acrt = np.empty(len(modelled)) + for j in range(len(modelled)): + if j == 0: + acrt[0] = att[0] / (modelled[0] - d_L) # backward diff to d_L + else: + acrt[j] = (att[j] - att[j - 1]) / (modelled[j] - modelled[j - 1]) + treated = udose - d_L > 1e-9 + overall = (dy[treated] - mu0).mean() + # Plug-in overall ACRT: density-weighted mean over treated units' doses. + tdose = udose[treated] + acrt_by_level = {d: acrt[j] for j, d in enumerate(modelled)} + overall_acrt = np.mean([acrt_by_level[min(modelled, key=lambda m: abs(m - d))] for d in tdose]) + n_c = int(cm.sum()) + se = np.empty(len(modelled)) + for j, d in enumerate(modelled): + tmask = np.abs(udose - d) <= 1e-9 + n_j = int(tmask.sum()) + it = (dy[tmask] - mu0 - att[j]) / n_j + ic = -(dy[cm] - mu0) / n_c + se[j] = np.sqrt((it**2).sum() + (ic**2).sum()) + return modelled, att, acrt, overall, overall_acrt, se + + class TestDiscreteSaturated: """Saturated regression for discrete/multi-valued treatment (CGBS 2024 Eq 4.1).""" @@ -1324,3 +1368,190 @@ def test_coverage(self, ci_params): assert np.all( (rate >= lo_band) & (rate <= 0.995) ), f"{'boot' if use_boot else 'analytic'} coverage {rate} off nominal" + + +class TestLowestDose: + """Remark 3.1 lowest-dose-as-control (CGBS 2024): estimand ATT(d) - ATT(d_L).""" + + _KW = dict( + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + dose="dose", + aggregate="dose", + ) + + def test_dL_to_zero_exact_equivalence(self): + """Relabelling the D=0 group as a tiny common dose d_L=eps reproduces the + never_treated fit EXACTLY (same control units, same dY) — for any eps.""" + effects = {1.0: 0.5, 2.0: 1.5, 4.0: 2.5} + df0 = _make_discrete_panel(effects, n_per_level=45, n_control=90, seed=11) + res_nt = ContinuousDiD(control_group="never_treated", treatment_type="discrete").fit( + df0, **self._KW + ) + eps = 1e-6 + df_re = df0.copy() + m = df_re["first_treat"] == 0 + df_re.loc[m, "first_treat"] = 1 + df_re.loc[m, "dose"] = eps + with warnings.catch_warnings(): + warnings.simplefilter("ignore") # tiny-boundary-gap warning + res_ld = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( + df_re, **self._KW + ) + # ATT and SE are EXACTLY equal (the differencing arithmetic is identical). + assert np.allclose( + res_ld.dose_response_att.effects, res_nt.dose_response_att.effects, atol=1e-10 + ) + assert np.allclose(res_ld.dose_response_att.se, res_nt.dose_response_att.se, atol=1e-10) + assert res_ld.reference_dose == eps + # Boundary ACRT references d_L=eps: ACRT(d_1) = ATT(d_1)/(d_1 - eps). + att1 = res_ld.dose_response_att.effects[0] + assert np.isclose(res_ld.dose_response_acrt.effects[0], att1 / (1.0 - eps), atol=1e-10) + + def test_hand_calc_att_acrt_overall(self): + """ATT(d)=mean(dY|d)-mean(dY|d_L); ACRT backward-to-d_L; overall = mean_T.""" + all_levels = [1.0, 2.0, 4.0] + df = _make_discrete_panel({1.0: 0.5, 2.0: 1.5, 4.0: 2.5}, n_control=0, seed=12) + res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( + df, **self._KW + ) + modelled, att, acrt, overall, overall_acrt, _ = _hand_calc_lowest_dose(df, all_levels) + assert np.allclose(res.dose_response_att.dose_grid, modelled) # d_L excluded + assert np.allclose(res.dose_response_att.effects, att, atol=1e-12) + assert np.allclose(res.dose_response_acrt.effects, acrt, atol=1e-12) + assert np.isclose(res.overall_att, overall, atol=1e-12) + # overall_acrt (plug-in scalar) asserted explicitly, not just the rows. + assert np.isclose(res.overall_acrt, overall_acrt, atol=1e-12) + + def test_hand_calc_analytical_se(self): + """Lowest-dose SE == per-level 2x2 DiD SE with the d_L group as control.""" + all_levels = [1.0, 2.0, 4.0] + df = _make_discrete_panel({1.0: 0.5, 2.0: 1.5, 4.0: 2.5}, n_control=0, seed=13) + res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( + df, **self._KW + ) + _, _, _, _, _, se = _hand_calc_lowest_dose(df, all_levels) + assert np.allclose(res.dose_response_att.se, se, atol=1e-10) + assert np.all(np.isfinite(res.dose_response_att.se)) + + def test_acrt_boundary_backward_to_dL(self): + """ACRT(d_1) = ATT(d_1)/(d_1 - d_L); ACRT(d_j>=2) adjacent backward diffs.""" + df = _make_discrete_panel({1.0: 0.5, 2.0: 1.5, 4.0: 2.5}, n_control=0, seed=14) + res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( + df, **self._KW + ) + att = res.dose_response_att.effects + acrt = res.dose_response_acrt.effects + modelled = res.dose_response_att.dose_grid # [2, 4]; d_L = 1 + assert np.isclose(acrt[0], att[0] / (modelled[0] - 1.0), atol=1e-12) + assert np.isclose(acrt[1], (att[1] - att[0]) / (modelled[1] - modelled[0]), atol=1e-12) + + def test_binary_above_reference(self): + """Single modelled dose above d_L (J=1): ACRT(d_1) = ATT(d_1)/(d_1 - d_L).""" + df = _make_discrete_panel({1.0: 0.4, 3.0: 2.0}, n_per_level=60, n_control=0, seed=15) + res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( + df, **self._KW + ) + assert res.reference_dose == 1.0 + assert np.allclose(res.dose_response_att.dose_grid, [3.0]) + att1 = res.dose_response_att.effects[0] + assert np.isclose(res.dose_response_acrt.effects[0], att1 / (3.0 - 1.0), atol=1e-12) + + def test_dgp_recovery_discrete(self): + """Recover ATT(d) - ATT(d_L) with known effects (no noise -> exact).""" + effects = {1.0: 0.5, 2.0: 2.0, 4.0: 1.0} # non-monotone + df = _make_discrete_panel(effects, n_per_level=60, n_control=0, noise=0.0, seed=16) + res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( + df, **self._KW + ) + # d_L = 1 (effect 0.5); ATT(2)=2.0-0.5=1.5, ATT(4)=1.0-0.5=0.5. + assert np.allclose(res.dose_response_att.effects, [1.5, 0.5], atol=1e-10) + # ACRT: bwd@2 = 1.5/(2-1) = 1.5; bwd@4 = (0.5-1.5)/(4-2) = -0.5. + assert np.allclose(res.dose_response_acrt.effects, [1.5, -0.5], atol=1e-10) + + def test_continuous_mass_point_recovery(self): + """Continuous B-spline with a mass point at d_L recovers a linear slope.""" + rng = np.random.default_rng(17) + rows, uid = [], 0 + beta = 0.7 # ATT(d) - ATT(d_L) = beta * (d - d_L) + for d in [1.0] * 70 + list(rng.uniform(1.5, 5.0, 210)): + base = rng.normal(0, 1) + for p in (0, 1): + y = base + 0.3 * p + (beta * (d - 1.0) if p >= 1 else 0.0) + rows.append((uid, p, y, 1, d)) + uid += 1 + df = pd.DataFrame(rows, columns=["unit", "period", "outcome", "first_treat", "dose"]) + res = ContinuousDiD(control_group="lowest_dose", treatment_type="continuous", degree=1).fit( + df, **self._KW + ) + assert res.reference_dose == 1.0 + assert np.all(res.dose_response_att.dose_grid > 1.0) # grid above d_L + grid = res.dose_response_att.dose_grid + expected = beta * (grid - 1.0) + # B-spline + finite sample: modest tolerance on the recovered curve. + assert np.allclose(res.dose_response_att.effects, expected, atol=0.15) + + def test_analytical_vs_bootstrap_se(self, ci_params): + """Analytical and multiplier-bootstrap SE agree (both carry d_L variance).""" + df = _make_discrete_panel({1.0: 0.5, 2.0: 1.5, 4.0: 2.5}, n_control=0, noise=1.0, seed=18) + res_a = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( + df, **self._KW + ) + nb = ci_params.bootstrap(999, min_n=199) + res_b = ContinuousDiD( + control_group="lowest_dose", treatment_type="discrete", n_bootstrap=nb, seed=7 + ).fit(df, **self._KW) + thr = 0.40 if nb < 200 else 0.15 + rel = ( + np.abs(res_b.dose_response_att.se - res_a.dose_response_att.se) + / res_a.dose_response_att.se + ) + assert np.all(rel < thr), f"boot/analytic SE disagree: {rel}" + + def test_pre_period_placebo(self): + """Pre-treatment cell att_glob ~ 0 (both treated and d_L control untreated).""" + effects = {1.0: 0.5, 2.0: 1.5, 4.0: 2.5} + df = _make_discrete_panel( + effects, n_per_level=50, n_control=0, noise=0.0, seed=19, cohorts=(2,), n_periods=3 + ) + res = ContinuousDiD( + control_group="lowest_dose", treatment_type="discrete", base_period="varying" + ).fit(df, **self._KW) + # Pre-period (t=1 < g=2) cell effects difference out to ~0. + pre = [v for (g, t), v in res.group_time_effects.items() if t < g] + assert pre, "expected a pre-period cell" + assert all(abs(c["att_glob"]) < 1e-9 for c in pre) + + @pytest.mark.slow + def test_coverage(self, ci_params): + """Analytical & bootstrap SE achieve nominal coverage for ATT(d)-ATT(d_L).""" + modelled = [2.0, 4.0] + effects = {1.0: 0.5, 2.0: 1.5, 4.0: 2.5} + truth = {d: effects[d] - effects[1.0] for d in modelled} # ATT(d) - ATT(d_L) + reps = 150 + for use_boot in (False, True): + nb = ci_params.bootstrap(299, min_n=99) if use_boot else 0 + lo_band = 0.85 if (use_boot and nb < 200) else 0.88 + cover = np.zeros(len(modelled)) + for s in range(reps): + df = _make_discrete_panel( + effects, n_per_level=45, n_control=0, noise=1.0, seed=40_000 + s + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = ContinuousDiD( + control_group="lowest_dose", + treatment_type="discrete", + n_bootstrap=nb, + seed=s, + ).fit(df, **self._KW) + lo = res.dose_response_att.conf_int_lower + hi = res.dose_response_att.conf_int_upper + for j, d in enumerate(modelled): + cover[j] += int(lo[j] <= truth[d] <= hi[j]) + rate = cover / reps + assert np.all( + (rate >= lo_band) & (rate <= 0.995) + ), f"{'boot' if use_boot else 'analytic'} coverage {rate} off nominal"