pyfector.diagnostics

  1"""
  2Diagnostic tests on a fitted :class:`~pyfector.fect.FectResult`.
  3
  4:func:`run_diagnostics` (also available as ``result.diagnose(...)``)
  5runs complementary checks on the pre-treatment portion of the event
  6study and on the post-treatment fit.
  7
  8* **Pre-trend F-test.** Wald-type joint significance of the
  9  pre-treatment ATTs against zero; small p-values cast doubt on the
 10  parallel-trends assumption.
 11* **Equivalence F-test.** Non-central F-test for whether all
 12  pre-treatment ATTs lie inside a user-supplied equivalence bound
 13  (``f_threshold``); small p-values reject inequivalence. Its noncentrality
 14  is ``N * f_threshold``, where ``N`` is the largest observed treated-unit
 15  count among the tested pre-periods. Bootstrap draws estimate covariance.
 16* **TOST** (two one-sided t-tests) on each individual pre-treatment
 17  coefficient with bound ``tost_threshold``. The default convention
 18  follows Liu et al. (2024): ``0.36 * sqrt(sigma2_fect)``, where
 19  ``sigma2_fect`` is the residual variance of the requested additive
 20  fixed-effect baseline (not of the chosen MC/IFE/CFE estimator). Pass
 21  an explicit positive float for an absolute outcome-scale bound.
 22* **Placebo test.** Select pre-treatment placebo periods, remove those
 23  observed cells from the model-fitting mask, refit with the selected
 24  model configuration, and estimate the placebo ATT from the withheld
 25  observations. This is the Liu-Wang-Xu/R-``fect`` holdout/refit
 26  placebo procedure, not an average of already-estimated pre-period
 27  event-study coefficients.
 28* **Carryover test.** Look for residual effect after treatment turns
 29  off (reversal units only). Equivalence p deferred until
 30  ``att_off_boot`` is plumbed through the inference layer.
 31* **Leave-one-out (LOO).** Drop each post-treatment period in turn and
 32  report the maximum change in overall ATT. This is post-treatment
 33  sensitivity, not the paper's leave-one-pre-period-out no-pretrend test.
 34
 35Result shape
 36------------
 37``Diagnostics`` is a slim-safe registry of per-test result objects.
 38Known tests remain available through ergonomic properties such as
 39``diag.tost`` and ``diag.pretrend_f``, while future tests can be added
 40under ``diag.tests[name]`` without changing the envelope schema. The
 41old ``DiagnosticResult`` name remains as a compatibility alias.
 42"""
 43
 44from __future__ import annotations
 45
 46from dataclasses import dataclass, field
 47from typing import Any
 48import warnings
 49
 50import numpy as np
 51
 52
 53@dataclass(frozen=True)
 54class TostResult:
 55    """Per-period TOST equivalence test on pre-treatment coefficients."""
 56    pvals: np.ndarray
 57    periods: np.ndarray
 58    threshold: float
 59    threshold_source: str
 60    sigma2_fect: float | None
 61    max_pval: float
 62    all_pass: bool
 63
 64    def summary(self) -> str:
 65        lines = [
 66            f"TOST per pre-period (threshold={self.threshold:.4f}, "
 67            f"source={self.threshold_source}):"
 68        ]
 69        for i, t in enumerate(self.periods):
 70            p = float(self.pvals[i])
 71            status = "PASS" if p < 0.05 else "fail"
 72            lines.append(f"  t={t:+.0f}: p={p:.4f} [{status}]")
 73        lines.append(
 74            f"  max p={self.max_pval:.4f}  all_pass={self.all_pass}"
 75        )
 76        return "\n".join(lines)
 77
 78
 79@dataclass(frozen=True)
 80class PretrendFResult:
 81    """Joint F-test for all pre-treatment ATTs equal to zero."""
 82    f_stat: float
 83    p_value: float
 84    df1: int
 85    df2: int
 86
 87    def summary(self) -> str:
 88        return (
 89            f"Pre-trend F-test:\n"
 90            f"  F({self.df1},{self.df2}) = {self.f_stat:.4f}, "
 91            f"p = {self.p_value:.4f}"
 92        )
 93
 94
 95@dataclass(frozen=True)
 96class EquivFResult:
 97    """Equivalence F-test (non-central F) for pre-trend bound."""
 98    p_value: float
 99    f_threshold: float
100    n_treated: int | None = None
101    ncp: float | None = None
102
103    def summary(self) -> str:
104        return (
105            f"Equivalence F-test:\n"
106            f"  p = {self.p_value:.4f} (f_threshold={self.f_threshold})"
107        )
108
109
110@dataclass(frozen=True)
111class PlaceboResult:
112    """Holdout/refit placebo ATT with bootstrap and TOST p-values.
113
114    ``n_obs`` counts withheld cells from eligible units. ``n_dropped_units``
115    counts units excluded after holdout; ``n_boot`` and ``n_boot_failed``
116    count valid and failed placebo bootstrap replications.
117    """
118    estimate: float
119    se: float
120    p_value: float
121    equiv_p_value: float
122    period: tuple[int, int]
123    n_obs: int
124    n_boot: int
125    n_dropped_units: int = 0
126    n_boot_failed: int = 0
127
128    def summary(self) -> str:
129        return (
130            f"Placebo test (holdout/refit, window {self.period}):\n"
131            f"  ATT={self.estimate:.4f} (SE={self.se:.4f}) "
132            f"p={self.p_value:.4f} equiv_p={self.equiv_p_value:.4f}"
133        )
134
135
136@dataclass(frozen=True)
137class CarryoverResult:
138    """Carryover-window mean ATT (post-reversal). SE/p deferred."""
139    estimate: float
140    period: tuple[int, int]
141
142    def summary(self) -> str:
143        return (
144            f"Carryover test (window {self.period}):\n"
145            f"  ATT={self.estimate:.4f}"
146        )
147
148
149@dataclass(frozen=True)
150class LooResult:
151    """Post-treatment sensitivity from dropping each post-treatment period.
152
153    This is not a leave-one-pre-period-out no-pretrend test.
154    """
155    atts: np.ndarray
156    periods: np.ndarray
157    max_change: float
158
159    def summary(self) -> str:
160        return (
161            f"Leave-one-out:\n"
162            f"  Max ATT change = {self.max_change:.6f}"
163        )
164
165
166_DIAGNOSTIC_SUMMARY_ORDER = (
167    "pretrend_f", "equiv_f", "tost", "placebo", "carryover", "loo",
168)
169
170
171@dataclass
172class Diagnostics:
173    """Slim-safe registry of diagnostic test results.
174
175    ``tests`` is the future-proof extension point. Built-in tests expose
176    convenience properties for stable user code, but the container itself
177    does not need a new field every time pyfector adds a diagnostic.
178    """
179    options: dict = field(default_factory=dict)
180    tests: dict[str, Any] = field(default_factory=dict)
181
182    @property
183    def available(self) -> tuple[str, ...]:
184        """Names of populated diagnostics, with built-ins first."""
185        known = [name for name in _DIAGNOSTIC_SUMMARY_ORDER if name in self.tests]
186        extra = sorted(name for name in self.tests if name not in set(known))
187        return tuple(known + extra)
188
189    def get(self, name: str, default: Any = None) -> Any:
190        """Return a diagnostic by name, or ``default`` if absent."""
191        return self.tests.get(name, default)
192
193    def set_test(self, name: str, value: Any | None) -> None:
194        """Set or remove a diagnostic result by name."""
195        if not name or not isinstance(name, str):
196            raise ValueError("diagnostic name must be a non-empty string")
197        if value is None:
198            self.tests.pop(name, None)
199        else:
200            self.tests[name] = value
201
202    def set(self, name: str, value: Any | None) -> None:
203        """Compatibility shortcut for :meth:`set_test`."""
204        self.set_test(name, value)
205
206    def __contains__(self, name: str) -> bool:
207        return name in self.tests
208
209    def __getitem__(self, name: str) -> Any:
210        return self.tests[name]
211
212    def _typed(self, name: str, typ: type) -> Any | None:
213        value = self.tests.get(name)
214        return value if isinstance(value, typ) else None
215
216    @property
217    def tost(self) -> TostResult | None:
218        return self._typed("tost", TostResult)
219
220    @tost.setter
221    def tost(self, value: TostResult | None) -> None:
222        self.set_test("tost", value)
223
224    @property
225    def pretrend_f(self) -> PretrendFResult | None:
226        return self._typed("pretrend_f", PretrendFResult)
227
228    @pretrend_f.setter
229    def pretrend_f(self, value: PretrendFResult | None) -> None:
230        self.set_test("pretrend_f", value)
231
232    @property
233    def equiv_f(self) -> EquivFResult | None:
234        return self._typed("equiv_f", EquivFResult)
235
236    @equiv_f.setter
237    def equiv_f(self, value: EquivFResult | None) -> None:
238        self.set_test("equiv_f", value)
239
240    @property
241    def placebo(self) -> PlaceboResult | None:
242        return self._typed("placebo", PlaceboResult)
243
244    @placebo.setter
245    def placebo(self, value: PlaceboResult | None) -> None:
246        self.set_test("placebo", value)
247
248    @property
249    def carryover(self) -> CarryoverResult | None:
250        return self._typed("carryover", CarryoverResult)
251
252    @carryover.setter
253    def carryover(self, value: CarryoverResult | None) -> None:
254        self.set_test("carryover", value)
255
256    @property
257    def loo(self) -> LooResult | None:
258        return self._typed("loo", LooResult)
259
260    @loo.setter
261    def loo(self, value: LooResult | None) -> None:
262        self.set_test("loo", value)
263
264    def summary(self) -> str:
265        lines = ["Diagnostic Tests", "=" * 50]
266        for name in self.available:
267            sub = self.tests[name]
268            if sub is not None:
269                if hasattr(sub, "summary"):
270                    lines.append(sub.summary())
271                else:
272                    lines.append(f"{name}: {sub!r}")
273        return "\n".join(lines)
274
275
276DiagnosticResult = Diagnostics
277
278
279_VALID_DIAG_NAMES = {"tost", "pretrend_f", "equiv_f", "placebo", "carryover", "loo"}
280_VALID_OPTION_KEYS = {
281    "tost_threshold", "f_threshold", "alpha",
282    "placebo_period", "carryover_period", "loo",
283}
284_REQUIRED_TESTS = ("tost", "pretrend_f", "equiv_f")
285
286
287def _is_positive_finite_number(value: Any) -> bool:
288    return (
289        isinstance(value, (int, float, np.integer, np.floating))
290        and not isinstance(value, (bool, np.bool_))
291        and np.isfinite(value)
292        and value > 0
293    )
294
295
296def _validate_tost_threshold_value(tost_threshold: Any) -> None:
297    if tost_threshold is None:
298        raise ValueError(
299            "tost_threshold=None is not a valid value. Omit the argument "
300            "for the Liu et al. (2024) default (0.36 * sqrt(sigma2_fect)) "
301            "or pass a positive float for an explicit absolute bound."
302        )
303    if not _is_positive_finite_number(tost_threshold):
304        raise ValueError(
305            f"tost_threshold must be a positive finite float, got "
306            f"{tost_threshold!r}"
307        )
308
309
310def _validate_period_option(name: str, value: Any) -> tuple[int, int]:
311    if (
312        not isinstance(value, tuple)
313        or len(value) != 2
314        or isinstance(value[0], bool)
315        or isinstance(value[1], bool)
316        or not np.isfinite(value[0])
317        or not np.isfinite(value[1])
318    ):
319        raise ValueError(
320            f"diagnostics_options['{name}'] must be a (start, end) tuple "
321            "of finite numeric periods."
322        )
323    start = int(value[0])
324    end = int(value[1])
325    if start > end:
326        raise ValueError(
327            f"diagnostics_options['{name}'] must have start <= end."
328        )
329    return start, end
330
331
332def _validate_placebo_period(value: Any) -> tuple[int, int]:
333    start, end = _validate_period_option("placebo_period", value)
334    if start >= 0:
335        raise ValueError(
336            "diagnostics_options['placebo_period'] must include at least "
337            "one pre-treatment relative period (< 0)."
338        )
339    if end > 0:
340        raise ValueError(
341            "diagnostics_options['placebo_period'] must not extend into "
342            "post-treatment periods (> 0)."
343        )
344    return start, end
345
346
347def _resolve_tost_threshold(result, tost_threshold) -> tuple[float, str, float | None]:
348    """Resolve the TOST equivalence bound.
349
350    Returns (resolved_absolute_value, source_label, sigma2_fect_used).
351    """
352    if tost_threshold == 0.36:
353        sigma2_fect = getattr(result, "sigma2_fect", 0.0) or 0.0
354        if sigma2_fect <= 0.0:
355            raise ValueError(
356                "result.sigma2_fect is missing or zero; refit with "
357                "pyfector >= 0.2.0 so sigma2_fect is populated, or pass "
358                "tost_threshold=<positive float> for an explicit absolute "
359                "bound."
360            )
361        return 0.36 * float(np.sqrt(sigma2_fect)), "auto", float(sigma2_fect)
362
363    _validate_tost_threshold_value(tost_threshold)
364
365    return float(tost_threshold), "explicit", None
366
367
368def run_diagnostics(
369    result,
370    f_threshold: float = 0.5,
371    tost_threshold: float = 0.36,
372    placebo_period: tuple[int, int] | None = None,
373    carryover_period: tuple[int, int] | None = None,
374    loo: bool = False,
375    alpha: float = 0.05,
376    *,
377    _requested: list[str] | None = None,
378    _strict: bool | None = None,
379) -> Diagnostics:
380    """Run diagnostic tests on a FectResult.
381
382    Parameters
383    ----------
384    result : FectResult
385        Must have inference results (``se=True``).
386    f_threshold : float
387        Noncentrality per treated unit for the equivalence F-test:
388        ``ncp = N * f_threshold``. The default coefficient is 0.5; in the
389        paper's notation this is ``kappa**2``. Small p-values support
390        equivalence. ``N`` is the largest ``count_on`` in the tested
391        pre-periods and also determines both F-tests' scale and degrees
392        of freedom.
393    tost_threshold : float
394        Equivalence bound for TOST. The literal ``0.36`` (default or
395        explicit) triggers Liu et al. (2024)'s scale-aware bound,
396        ``0.36 * sqrt(result.sigma2_fect)``. Any other positive float
397        is taken as an absolute outcome-scale bound. ``None`` is
398        invalid.
399    placebo_period : (start, end), optional
400        Relative pre-treatment window for the holdout/refit placebo test.
401        Observed cells with ``start <= time_on <= end`` and ``time_on < 0``
402        are removed from the fitting mask, the model is refit using the
403        selected rank/lambda configuration, and the placebo ATT is
404        computed from those withheld cells. For example, ``(-3, 0)``
405        withholds relative periods ``-3, -2, -1``.
406        Treated and reversal units must retain at least ``max(min_T0, 1)``
407        untreated observations after holdout. Bootstrap draws use the
408        fit's cluster column when supplied.
409    carryover_period : (start, end), optional
410        Relative time window for carryover test.
411    loo : bool
412        If True, run leave-one-out post-period sensitivity.
413    alpha : float
414        Significance cutoff used for ``TostResult.all_pass``.
415
416    Returns
417    -------
418    Diagnostics
419        Available tests and their resolved options. Unavailable configured
420        tests raise for an explicit list and warn in best-effort mode.
421    """
422    from scipy import stats
423
424    diag = Diagnostics()
425
426    requested = set(_requested) if _requested is not None else _VALID_DIAG_NAMES.copy()
427    strict = _requested is not None if _strict is None else _strict
428    expected = requested.copy()
429    unavailable_reasons = {}
430    if _requested is None:
431        if placebo_period is None:
432            expected.discard("placebo")
433        if carryover_period is None:
434            expected.discard("carryover")
435        if not loo:
436            expected.discard("loo")
437
438    def report_unavailable():
439        missing = sorted((expected - set(diag.tests)) | (expected & set(unavailable_reasons)))
440        if missing:
441            descriptions = [
442                f"{name} ({unavailable_reasons[name]})" if name in unavailable_reasons else name
443                for name in missing
444            ]
445            message = f"Requested diagnostics unavailable: {', '.join(descriptions)}"
446            if strict:
447                raise ValueError(message)
448            warnings.warn(message, RuntimeWarning, stacklevel=3)
449
450    if result.inference is None:
451        report_unavailable()
452        return diag
453
454    inf = result.inference
455    time_on = result.time_on
456    att_on = result.att_on
457
458    if time_on is None or att_on is None:
459        report_unavailable()
460        return diag
461
462    pre_mask = time_on < 0
463
464    pre_idx = np.where(pre_mask)[0]
465    k = len(pre_idx)
466    att_pre = att_on[pre_idx]
467
468    if placebo_period is not None:
469        placebo_period = _validate_placebo_period(placebo_period)
470
471    needs_tost_threshold = bool({"tost", "placebo"} & requested)
472    threshold_abs = float("nan")
473    threshold_source: str | None = None
474    sigma2_fect_used: float | None = None
475    if needs_tost_threshold:
476        _validate_tost_threshold_value(tost_threshold)
477        try:
478            threshold_abs, threshold_source, sigma2_fect_used = _resolve_tost_threshold(
479                result, tost_threshold
480            )
481        except ValueError as exc:
482            if strict or _strict is None:
483                raise
484            for name in requested & {"tost", "placebo"}:
485                unavailable_reasons[name] = str(exc)
486
487    # Pre-trend / equivalence F-tests share the bootstrap covariance pass.
488    if (
489        k > 0 and ("pretrend_f" in requested or "equiv_f" in requested)
490        and inf.att_on_boot is not None
491    ):
492        boot_pre_all = inf.att_on_boot[pre_idx, :]
493        valid_boot = np.all(np.isfinite(boot_pre_all), axis=0)
494        boot_pre = boot_pre_all[:, valid_boot]
495        n_boot = boot_pre.shape[1]
496        if n_boot <= k:
497            boot_pre = None
498
499    else:
500        boot_pre = None
501
502    if boot_pre is not None:
503        counts = getattr(result, "count_on", None)
504        if counts is None:
505            reason = "observed treated-unit counts (count_on) are required"
506            n_treated = 0
507        else:
508            pre_counts = np.asarray(counts)[pre_idx]
509            valid_counts = np.all(np.isfinite(pre_counts)) and np.all(pre_counts >= 0)
510            n_treated = int(np.max(pre_counts)) if valid_counts else 0
511            reason = "insufficient observed treated units for the number of tested pre-periods"
512        if n_treated <= k:
513            for name in requested & {"pretrend_f", "equiv_f"}:
514                unavailable_reasons[name] = reason
515            boot_pre = None
516
517    if boot_pre is not None:
518        S = np.cov(boot_pre)
519        if k == 1:
520            S = np.asarray(S).reshape(1, 1)
521
522        try:
523            cond = np.linalg.cond(S)
524            S_inv = np.linalg.pinv(S) if cond > 1e12 else np.linalg.inv(S)
525            F_raw = float(att_pre @ S_inv @ att_pre)
526            scale = (n_treated - k) / ((n_treated - 1) * k)
527            F_stat = F_raw * scale
528
529            if np.isfinite(F_stat) and F_stat >= 0:
530                if "pretrend_f" in requested:
531                    diag.pretrend_f = PretrendFResult(
532                        f_stat=F_stat,
533                        p_value=float(stats.f.sf(F_stat, k, n_treated - k)),
534                        df1=k,
535                        df2=n_treated - k,
536                    )
537                if "equiv_f" in requested:
538                    ncp = n_treated * f_threshold
539                    diag.equiv_f = EquivFResult(
540                        p_value=float(stats.ncf.cdf(F_stat, k, n_treated - k, ncp)),
541                        f_threshold=f_threshold,
542                        n_treated=n_treated,
543                        ncp=ncp,
544                    )
545        except np.linalg.LinAlgError:
546            pass
547
548    # Per-period TOST.
549    if k > 0 and "tost" in requested and "tost" not in unavailable_reasons and inf.att_on_se is not None:
550        se_pre = inf.att_on_se[pre_idx]
551
552        tost_pvals = np.full(k, np.nan)
553        for i in range(k):
554            if se_pre[i] > 0:
555                if inf.att_on_boot is not None:
556                    n_boot_i = int(np.isfinite(inf.att_on_boot[pre_idx[i], :]).sum())
557                else:
558                    n_boot_i = 200
559                df = max(n_boot_i - 1, 1)
560                t_upper = (att_pre[i] - threshold_abs) / se_pre[i]
561                t_lower = (att_pre[i] + threshold_abs) / se_pre[i]
562                p_upper = float(stats.t.cdf(t_upper, df))
563                p_lower = float(1 - stats.t.cdf(t_lower, df))
564                tost_pvals[i] = max(p_upper, p_lower)
565
566        finite = tost_pvals[np.isfinite(tost_pvals)]
567        max_p = float(finite.max()) if finite.size else float("nan")
568        all_pass = bool(finite.size == k and (finite < alpha).all())
569        if finite.size < k:
570            unavailable_reasons["tost"] = "missing usable standard errors for pre-treatment periods"
571        if finite.size:
572            diag.tost = TostResult(
573                pvals=tost_pvals,
574                periods=time_on[pre_idx].copy(),
575                threshold=threshold_abs,
576                threshold_source=threshold_source,
577                sigma2_fect=sigma2_fect_used,
578                max_pval=max_p,
579                all_pass=all_pass,
580            )
581
582    # Holdout/refit placebo test. This follows R fect's placeboTest path:
583    # selected pre-treatment cells are removed from II before fitting.
584    if "placebo" in requested and "placebo" not in unavailable_reasons and placebo_period is not None:
585        try:
586            diag.placebo = _run_placebo_test(
587                result,
588                placebo_period=placebo_period,
589                threshold_abs=threshold_abs,
590            )
591        except (ValueError, RuntimeError) as exc:
592            if strict or _strict is None:
593                raise
594            unavailable_reasons["placebo"] = str(exc)
595
596    # Carryover (estimate only; bootstrap p deferred).
597    if (
598        "carryover" in requested
599        and carryover_period is not None
600        and hasattr(result, "att_off")
601        and result.att_off is not None
602    ):
603        c_start, c_end = carryover_period
604        time_off = getattr(result, "time_off", None)
605        if time_off is not None:
606            off_mask = (time_off >= c_start) & (time_off <= c_end)
607            if np.any(off_mask):
608                diag.carryover = CarryoverResult(
609                    estimate=float(np.mean(result.att_off[off_mask])),
610                    period=(int(c_start), int(c_end)),
611                )
612
613    # Leave-one-out post-period.
614    if "loo" in requested and loo:
615        post_mask = time_on >= 0
616        post_idx = np.where(post_mask)[0]
617        if len(post_idx) > 1:
618            full_att = result.att_avg
619            loo_atts: list[float] = []
620            loo_periods: list[float] = []
621            for drop_i in post_idx:
622                remaining = np.delete(post_idx, np.where(post_idx == drop_i))
623                if len(remaining) > 0:
624                    loo_atts.append(float(np.mean(att_on[remaining])))
625                    loo_periods.append(float(time_on[drop_i]))
626            atts_arr = np.array(loo_atts)
627            diag.loo = LooResult(
628                atts=atts_arr,
629                periods=np.array(loo_periods),
630                max_change=float(np.max(np.abs(atts_arr - full_att))),
631            )
632
633    diag.options = {
634        "tost_threshold": threshold_abs if needs_tost_threshold else None,
635        "tost_threshold_source": threshold_source,
636        "f_threshold": f_threshold,
637        "alpha": alpha,
638        "placebo_period": placebo_period,
639        "carryover_period": carryover_period,
640        "loo": loo,
641        "sigma2_fect": sigma2_fect_used if sigma2_fect_used is not None
642                       else getattr(result, "sigma2_fect", None),
643    }
644    try:
645        from . import __version__ as _ver
646        diag.options["pyfector_version"] = _ver
647    except Exception:
648        pass
649
650    report_unavailable()
651    return diag
652
653
654def _selected_boot_count(result) -> int:
655    inf = getattr(result, "inference", None)
656    boot = getattr(inf, "att_on_boot", None)
657    if boot is not None and getattr(boot, "ndim", 0) == 2 and boot.shape[1] > 0:
658        return int(boot.shape[1])
659    avg_boot = getattr(inf, "att_avg_boot", None)
660    if avg_boot is not None and len(avg_boot) > 0:
661        return int(len(avg_boot))
662    opts = getattr(result, "fit_options", {}) or {}
663    return int(opts.get("nboots", 200))
664
665
666def _estimate_with_selected_config(
667    method: str,
668    Y,
669    Y0,
670    X,
671    II,
672    W,
673    beta0,
674    *,
675    r_cv: int | None,
676    lambda_cv: float | None,
677    force_int: int,
678    tol: float,
679    max_iter: int,
680):
681    from .estimators import estimate_cfe, estimate_ife, estimate_mc
682
683    if method in ("fe", "ife", "both"):
684        return estimate_ife(
685            Y, Y0, X, II, W, beta0,
686            r=r_cv or 0, force=force_int, tol=tol, max_iter=max_iter,
687        )
688    if method == "mc":
689        return estimate_mc(
690            Y, Y0, X, II, W, beta0,
691            lam=0.0 if lambda_cv is None else lambda_cv,
692            force=force_int, tol=tol, max_iter=max_iter,
693        )
694    if method == "cfe":
695        return estimate_cfe(
696            Y, Y0, X, II, W, beta0,
697            r=r_cv or 0, force=force_int, tol=tol, max_iter=max_iter,
698        )
699    raise ValueError(f"Unknown method for placebo test: {method!r}")
700
701
702def _run_placebo_test(
703    result,
704    *,
705    placebo_period: tuple[int, int],
706    threshold_abs: float,
707) -> PlaceboResult:
708    """Run the holdout/refit placebo test used by Liu-Wang-Xu/R fect."""
709    from joblib import Parallel, delayed, parallel_config
710    from scipy import stats
711
712    from .backend import to_device, to_numpy
713    from .panel import initial_fit
714    from .inference import _make_bootstrap_sampler, _check_bootstrap_replicates
715
716    panel = getattr(result, "panel", None)
717    if panel is None:
718        raise ValueError(
719            "placebo_period requires result.panel because the placebo "
720            "test refits on a modified fitting mask."
721        )
722
723    opts = getattr(result, "fit_options", {}) or {}
724    if "force_int" not in opts:
725        raise ValueError(
726            "placebo_period requires fit metadata from pyfector.fect(). "
727            "Refit the model before running the placebo test."
728        )
729
730    p_start, p_end = placebo_period
731    holdout_mask = (
732        (panel.I > 0)
733        & np.isfinite(panel.T_on)
734        & (panel.T_on >= p_start)
735        & (panel.T_on <= p_end)
736        & (panel.T_on < 0)
737    )
738    if not np.any(holdout_mask):
739        raise ValueError(
740            "placebo_period selects no observed pre-treatment cells."
741        )
742
743    II_holdout = panel.II.copy()
744    II_holdout[holdout_mask] = 0.0
745
746    minimum = max(int(opts.get("min_T0", 1)), 1)
747    treated = panel.unit_type != 1
748    remaining = (II_holdout > 0).sum(axis=0)
749    eligible = ~treated | (remaining >= minimum)
750    if opts.get("min_T0_strict", False):
751        eligible &= remaining >= minimum
752    if not np.any(eligible & treated):
753        raise ValueError("No treated units remain eligible for the placebo refit after holdout")
754    n_dropped_units = int(np.count_nonzero(~eligible))
755    eligible_idx = np.flatnonzero(eligible)
756    holdout_mask[:, ~eligible] = False
757    n_obs = int(holdout_mask.sum())
758    if n_obs == 0:
759        raise ValueError("placebo_period selects no eligible observed pre-treatment cells")
760
761    normalize = bool(opts.get("normalize", False))
762    norm_factor = float(opts.get("norm_factor", 1.0) or 1.0)
763    if normalize and norm_factor != 1.0:
764        Y_np = panel.Y / norm_factor
765    else:
766        Y_np = panel.Y
767
768    force_int = int(opts["force_int"])
769    tol = float(opts.get("tol", 1e-7))
770    boot_tol = max(tol, 1e-3)
771    max_iter = int(opts.get("max_iter", 5000))
772    method = getattr(result, "method")
773    r_cv = getattr(result, "r_cv", None)
774    lambda_cv = getattr(result, "lambda_cv", None)
775
776    def _fit_placebo(unit_idx, tol_value):
777        Y_sub_np = Y_np[:, unit_idx]
778        II_sub_np = II_holdout[:, unit_idx]
779        X_sub_np = panel.X[:, unit_idx, :] if panel.X is not None else None
780        W_sub_np = panel.W[:, unit_idx] if panel.W is not None else None
781        mask_sub = holdout_mask[:, unit_idx]
782        if not np.any(mask_sub):
783            return None
784
785        Y_sub = to_device(Y_sub_np)
786        II_sub = to_device(II_sub_np)
787        X_sub = to_device(X_sub_np) if X_sub_np is not None else None
788        W_sub = to_device(W_sub_np) if W_sub_np is not None else None
789
790        Y0_sub, beta0_sub = initial_fit(Y_sub, X_sub, II_sub, force_int)
791        est = _estimate_with_selected_config(
792            method,
793            Y_sub,
794            Y0_sub,
795            X_sub,
796            II_sub,
797            W_sub,
798            beta0_sub,
799            r_cv=r_cv,
800            lambda_cv=lambda_cv,
801            force_int=force_int,
802            tol=tol_value,
803            max_iter=max_iter,
804        )
805        eff_sub = to_numpy(Y_sub - est.fit)
806        if normalize and norm_factor != 1.0:
807            eff_sub = eff_sub * norm_factor
808        return float(np.mean(eff_sub[mask_sub]))
809
810    point = _fit_placebo(eligible_idx, tol)
811    if point is None:
812        raise ValueError(
813            "placebo_period selects no observed pre-treatment cells."
814        )
815
816    nboots = _selected_boot_count(result)
817    n_jobs = int(opts.get("n_jobs", 1) or 1)
818    if n_jobs == 0:
819        n_jobs = 1
820
821    seed = getattr(result, "seed", None)
822    rng = np.random.default_rng(seed)
823    cluster = getattr(panel, "cluster", None)
824    sample_units = _make_bootstrap_sampler(
825        panel.unit_type[eligible_idx],
826        cluster[eligible_idx] if cluster is not None else None,
827    )
828
829    def _one_boot(boot_seed):
830        boot_rng = np.random.default_rng(boot_seed)
831        unit_idx = eligible_idx[sample_units(boot_rng)]
832        try:
833            return _fit_placebo(unit_idx, boot_tol)
834        except Exception:
835            return None
836
837    boot_seeds = rng.integers(0, 2**31, size=nboots)
838    if n_jobs == 1:
839        boot_results = [_one_boot(s) for s in boot_seeds]
840    else:
841        with parallel_config(backend="loky", inner_max_num_threads=1):
842            boot_results = Parallel(n_jobs=n_jobs)(
843                delayed(_one_boot)(s) for s in boot_seeds
844            )
845    boot = np.array(
846        [b for b in boot_results if b is not None and np.isfinite(b)],
847        dtype=float,
848    )
849    n_failed = _check_bootstrap_replicates(boot.size, nboots, "placebo bootstrap")
850
851    se = float("nan")
852    p_value = float("nan")
853    equiv_p = float("nan")
854    if boot.size > 1:
855        se = float(np.std(boot, ddof=1))
856        if se > 0:
857            # Match R fect's default placebo-test summary: a normal
858            # approximation from the bootstrap SE, not a percentile p-value.
859            p_value = float(2 * stats.norm.sf(abs(point / se)))
860            z_upper = (point - threshold_abs) / se
861            z_lower = (point + threshold_abs) / se
862            equiv_p = max(
863                float(stats.norm.cdf(z_upper)),
864                float(1 - stats.norm.cdf(z_lower)),
865            )
866
867    return PlaceboResult(
868        estimate=float(point),
869        se=se,
870        p_value=p_value,
871        equiv_p_value=equiv_p,
872        period=(int(p_start), int(p_end)),
873        n_obs=n_obs,
874        n_boot=int(boot.size),
875        n_dropped_units=n_dropped_units,
876        n_boot_failed=n_failed,
877    )
878
879
880def validate_diagnostics_request(
881    diagnostics: Any,
882    diagnostics_options: dict | None,
883    se: bool,
884) -> list[str] | None:
885    """Validate a fect()-side diagnostics request and return the resolved
886    list of test names to run (or ``None`` if ``diagnostics='none'``).
887
888    Raises ``ValueError`` for any invalid combination, all before any
889    expensive estimation runs.
890    """
891    opts = diagnostics_options or {}
892    bad_keys = set(opts.keys()) - _VALID_OPTION_KEYS
893    if bad_keys:
894        raise ValueError(
895            f"unknown diagnostics_options key(s): {sorted(bad_keys)}; "
896            f"valid: {sorted(_VALID_OPTION_KEYS)}"
897        )
898
899    if diagnostics == "none":
900        return None
901
902    if not se:
903        raise ValueError(
904            "diagnostics requires se=True; pre-trend F-test, TOST, and "
905            "placebo p-values all need bootstrap covariance."
906        )
907
908    if diagnostics == "full":
909        requested = list(_REQUIRED_TESTS)
910        if "placebo_period" in opts:
911            requested.append("placebo")
912        if "carryover_period" in opts:
913            requested.append("carryover")
914        if opts.get("loo", False):
915            requested.append("loo")
916    elif isinstance(diagnostics, list):
917        if len(diagnostics) == 0:
918            raise ValueError(
919                "diagnostics=[] is invalid; pass 'none' or a non-empty list."
920            )
921        bad = [d for d in diagnostics if d not in _VALID_DIAG_NAMES]
922        if bad:
923            raise ValueError(
924                f"Unknown diagnostic(s): {bad}; valid: "
925                f"{sorted(_VALID_DIAG_NAMES)}"
926            )
927        if "placebo" in diagnostics and "placebo_period" not in opts:
928            raise ValueError(
929                "diagnostics list includes 'placebo' but "
930                "diagnostics_options['placebo_period'] is missing."
931            )
932        if "carryover" in diagnostics and "carryover_period" not in opts:
933            raise ValueError(
934                "diagnostics list includes 'carryover' but "
935                "diagnostics_options['carryover_period'] is missing."
936            )
937        requested = list(diagnostics)
938    else:
939        raise ValueError(
940            f"diagnostics must be 'none', 'full', or a list of names; got "
941            f"{diagnostics!r}"
942        )
943
944    if "alpha" in opts:
945        alpha = opts["alpha"]
946        if (
947            not isinstance(alpha, (int, float, np.integer, np.floating))
948            or isinstance(alpha, (bool, np.bool_))
949            or not np.isfinite(alpha)
950            or not 0 < alpha < 1
951        ):
952            raise ValueError(
953                "diagnostics_options['alpha'] must be a finite float in (0, 1)."
954            )
955
956    if "f_threshold" in opts and not _is_positive_finite_number(opts["f_threshold"]):
957        raise ValueError(
958            "diagnostics_options['f_threshold'] must be a positive finite float."
959        )
960
961    if "tost_threshold" in opts and {"tost", "placebo"} & set(requested):
962        _validate_tost_threshold_value(opts["tost_threshold"])
963
964    if "placebo_period" in opts:
965        _validate_placebo_period(opts["placebo_period"])
966    if "carryover_period" in opts:
967        _validate_period_option("carryover_period", opts["carryover_period"])
968
969    if "loo" in opts and not isinstance(opts["loo"], (bool, np.bool_)):
970        raise ValueError("diagnostics_options['loo'] must be a bool.")
971
972    return requested