r"""Omitted-variable-bias sensitivity analysis.
Positivity diagnostics tell you whether the data can support the estimate. This
module answers a different question: *how strong would an unmeasured confounder
have to be to overturn the conclusion?*
Chernozhukov, Cinelli, Newey, Sharma & Syrgkanis (2026) bound the bias from omitting a
confounder by a product of three interpretable pieces. Their Theorem 2 gives the bias of a
linear functional, and Equation (14) in their Section 4 gives the bounds:
.. math::
|\mathrm{bias}| \le |\rho| \sqrt{\frac{c_D^2}{1 - c_D^2}}\; c_Y\;
\underbrace{\sqrt{\sigma^2 \nu^2}}_{\text{bias scale}}
where, writing :math:`\alpha` for the Riesz representer of the target parameter,
* :math:`\sigma^2 = E[(Y - \bar Q(A, W))^2]` -- residual outcome variance,
* :math:`\nu^2 = E[\alpha(A, W)^2]` -- how hard the estimand has to work to
extrapolate (it is large exactly when overlap is poor),
* ``cf_y`` -- the share of *residual* outcome variance the confounder would explain,
* ``cf_d`` -- the share of the Riesz representer's second moment that only the
confounder adds; for the ATE, the gain in precision of the treatment model,
* ``rho`` -- how adversarially aligned those two are; ``rho = 1`` is the worst case.
:math:`\sqrt{\sigma^2 \nu^2}` is a scale, not the largest bias: the multiplier
``|rho| sqrt(cf_y cf_d / (1 - cf_d))`` exceeds 1 for strong confounders (about 2.85 at
``cf_y = cf_d = 0.9``). The field that holds it keeps DoubleML's name, ``max_bias``.
Two things make this useful rather than merely formal. First, the *robustness
value* :func:`robustness_value` inverts the bound: it reports the single number
``cf_y = cf_d = RV`` at which the conclusion would flip, so there is no need to
guess sensitivity parameters at all. Second, :func:`benchmark` calibrates
``cf_y``/``cf_d`` against covariates you *did* measure -- "a confounder as strong as
age" is a claim a reader can evaluate, where "cf_y = 0.03" is not.
The parameterisation, including the definition of the benchmark gain statistics,
follows DoubleML's ``sensitivity_analysis``, so a strength means the same thing in both
libraries. The confidence limits of an ATT or an ATC differ: this module adds the
influence term of the conditioning share to the curve of ``nu^2``, and
``DoubleMLIRM._sensitivity_element_est`` omits it (doubleml-for-py at ``b69f86ef``). The
estimates differ too, because DoubleML fits a cross-fitted AIPW score and this package a
TMLE.
Inference. Each bound has one-sided limits, from Theorem 4 in Section 4 of Chernozhukov et
al. For the ATT and the ATC the curve of the doubly robust ``nu^2`` also carries
``-2 nu^2 (1{A = c} - p) / p``, the influence of the estimated share ``p`` of the
conditioning arm ``c``. The paper gives this kind of term for the estimate only (Online
Appendix A, "Statistical Inference", Equation (15)). :func:`_conditioning_share_influence`
holds this package's derivation.
The plug-in estimator of ``nu^2`` has no derived curve, so its bounds publish no limits
(docs/roadmap.md F26).
Scope: the bound applies to the linear functionals this library estimates -- the
counterfactual means, their contrasts against the reference arm, and the two conditional
effects. With more than two arms those are named for the arms they are about, so the
estimand to ask for is ``"ate[medium vs low]"`` rather than ``"ate"``; it is one bound
per contrast, because :math:`\nu^2` is the second moment of *that contrast's* Riesz
representer. Ratios are not linear functionals of the outcome regression, so use
:mod:`cleverly.sensitivity.evalue` for those.
Scope of the refusals: :data:`_FIT_WIDE_BOUND_RULES` refuses the fits below, for five
reasons.
* A DR-TMLE fit and a collaborative TMLE fit break a premise. The ``nu^2`` score of
Lemma 3 of Chernozhukov et al. and the standard error of their Theorem 4 assume a
consistently estimated treatment mechanism.
* A fit with a response mechanism, a fit with an intermediate variable, and a fit indexed
by a regime, a shift, or an MSM term are well posed. Theorem 2 of Chernozhukov et al.
covers each one as a linear functional of the outcome regression, and no implementation
is registered.
* An ``ipsi`` fit is outside that theorem, because its estimand contains the treatment
mechanism.
* A longitudinal fit has no registered derivation, and a median-combined repeated fit has
no influence function for the median bound.
* An arm-indexed fit that reports no counterfactual mean and no linear contrast has no
parameter to bound. A ratio-only fit, a population attributable risk or fraction, and
a complete-outcome natural-course mean are such fits.
Every entry point in this module reaches those rules through :func:`sensitivity_elements`,
and :class:`~cleverly.assessment.SensitivityFacade` declares the same reason on the
matching capability rows, so a fit the bound refuses is never advertised as available.
"""
from __future__ import annotations
import functools
from collections.abc import Callable
from dataclasses import dataclass, replace
from typing import TYPE_CHECKING, Any, ParamSpec
import numpy as np
from scipy import optimize
from .._inference_status import InferenceStatus, is_labelled
from .._refusal import NO_ALTERNATIVE, Refusal, RefusalReason
from .._support import record_refusal_site
from .._typing import FloatArray
from ..assessment import SENSITIVITY_ROUTES
from ..estimators.ctmle import CTMLEDerivativeFit, CTMLESelection, declares_no_selection
from ..estimators.direct_effect import declares_intermediate
from ..estimators.targeting import build_submodel
from ..exceptions import CapabilityError, DataError, refuse_inference, repeats_refusal
from ..inference.cluster import influence_variance
from ..inference.delta import one_sided_quantile
from ..inference.influence import spread_name
from ..interventions.learned import LEARNED_RULE_BOUND_REFUSAL
from ..targets import parameter_stem
from ..targets.population_intervention import natural_course_names, reads_natural_course_mean
from ..utils.bounds import g_bounds_for
from ..utils.random import resolve_assessment_seed
from ..utils.text import format_table
from ._parameters import ArmParameter, reported_arm_parameters, stratum_refusal
from .missingness import fit_wide_tilt_refusal
if TYPE_CHECKING: # pragma: no cover - typing only
from ..estimators._nuisance import RepeatFit
from ..estimators.base import TMLEResult
__all__ = [
"BenchmarkResult",
"SensitivityBounds",
"SensitivityElements",
"benchmark",
"omitted_variable_bounds",
"robustness_value",
"sensitivity_elements",
]
#: Targets for which the Riesz representer -- and therefore this bound -- is defined.
#:
#: Stems rather than reported names, since a multi-valued treatment names each parameter
#: for the arms it is about: ``ate[medium vs low]`` has a stem of ``ate`` and a
#: representer of its own. :func:`~cleverly.sensitivity._parameters.arm_parameters` is
#: what turns a fit's arms into the names these stems produce.
LINEAR_ESTIMANDS: frozenset[str] = frozenset({"ate", "ey", "ey1", "ey0", "att", "atc"})
#: The two ratio scales an E-value answers for, and the pointer a refusal appends when a
#: request or a fit names one. Offered for anything else, such as an attributable
#: fraction, the pointer sends the reader from one refusal to a second one.
_EVALUE_STEMS: frozenset[str] = frozenset({"rr", "or"})
_EVALUE_POINTER = " For a risk ratio or odds ratio use sensitivity.evalue()."
#: Accepted values of ``nu2_estimator``, in one place so the five docstrings that list
#: them and the refusal that rejects the rest cannot drift apart.
#: ``tests/unit/test_omitted_variable_refusals.py`` parses the ``{...}`` literal out of
#: every one of those docstrings and compares it with this tuple.
NU2_ESTIMATORS: tuple[str, ...] = ("auto", "doubly_robust", "plugin")
#: The capability rows :class:`~cleverly.assessment.SensitivityFacade` fills from this
#: module's refusal, named once so the tests that check those rows read the set rather
#: than respell it. ``benchmark`` shares the refusal on every point fit. On a
#: longitudinal fit it names its own missing derivation instead. Read off
#: :data:`~cleverly.assessment.SENSITIVITY_ROUTES` rather than listed here, because that
#: table already says which operations this module answers, and a second list is the
#: registry its comment warns against.
OMITTED_VARIABLE_OPERATIONS: tuple[str, ...] = tuple(
operation
for operation, route in SENSITIVITY_ROUTES.items()
if route.module == "omitted_variable"
)
#: What each non-arm parameter axis reports, named so the refusal says which functional
#: the bound is missing rather than which one the fit lacks. ``ipsi`` is absent on
#: purpose: an incremental intervention is refused for a different reason, and
#: :func:`_refuse_non_arm_axis` states it separately.
_AXIS_NOUNS: dict[str, str] = {
"regime": "A regime mean",
"policy": "A modified-policy mean",
"rr_tilt": "A risk-ratio tilt mean",
"msm": "A point-treatment MSM coefficient",
}
#: The longitudinal stop. ``tests/unit/test_assessment_contract.py`` matches the declared
#: row against the raised error, and both read this one record. The citation names the
#: research-queue row that holds the missing law.
_LONGITUDINAL_BOUND_REFUSAL = Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"no longitudinal sensitivity derivation is registered: the bound needs a per-period "
"latent law that identifies the bias of a sequential g-formula target "
"(docs/roadmap.md F16 tracks this stop).",
alternative=NO_ALTERNATIVE,
roadmap="F16",
)
#: The clause both mechanism refusals end on. Neither estimator assumes a consistent
#: treatment mechanism, and that assumption is what the published bound's ``nu^2`` and
#: its standard error are derived under.
_NO_NU2_DERIVATION = (
"No derivation registered here gives nu^2, or the bound's standard error, for an "
"estimator that does not assume a consistent treatment mechanism."
)
_DRTMLE_BOUND_REFUSAL = Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"the omitted-variable bound has no nu^2 estimate for a 'drtmle' fit. The default "
"estimator E[2 m(alpha_hat) - alpha_hat^2] equals nu_0^2 minus the squared error of "
"the fitted representer, by the Riesz identity, so it falls exactly where the fitted "
"mechanism is wrong, which is the case DR-TMLE guards against. The bound would be too "
"narrow and the robustness value too large. " + _NO_NU2_DERIVATION,
alternative="Fit the ordinary TMLE on the same data for the bound.",
roadmap="X50",
)
_CTMLE_BOUND_REFUSAL = Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"the omitted-variable bound has no nu^2 estimate for a 'collaborative_tmle' fit. The "
"working mechanism conditions on a function V of W: the selected adjustment set W_S "
"on the greedy, ordered and discrete paths, the selected candidate mechanism on the "
"lasso and learners paths, and the fitted outcome regression under "
"'oat'. Where the working mechanism is P(A | V) in the limit, the representer is "
"E[alpha_W | A, V], whose second moment cannot exceed the second moment of alpha_W, "
"while sigma^2 still comes from a regression on every "
"declared covariate. The product sigma^2 nu^2 belongs to no single conditioning set, "
"and the collaborative robustness value is optimistic by construction. " + _NO_NU2_DERIVATION,
alternative="Fit the ordinary TMLE on the same data for the bound.",
roadmap="X51",
)
#: The same refusal for the two collaborative fits that select nothing. Their mechanism
#: is fitted on every covariate, so the argument above does not apply to them, and the
#: sentence names the scope rule instead.
_CTMLE_FULL_CANDIDATE_BOUND_REFUSAL = Refusal(
RefusalReason.API_REDIRECTION,
"the omitted-variable bound is not offered on a 'collaborative_tmle' fit. This fit "
"selects nothing: a strategy='discrete' fit that declares one candidate, equal to the "
"full adjustment set, or a strategy='learners' fit that declares no candidate, equals "
"the ordinary TMLE fit. The bound reads the fitted method and has no collaborative "
"route.",
alternative="Fit the ordinary TMLE (TMLE, or TMLEMethod) on the same data for the bound.",
)
#: The clause the two well-posed mechanism refusals share: each functional is linear in
#: an outcome regression, which is the hypothesis of the general bound.
_THEOREM_2_COVERS = (
"so Theorem 2 of Chernozhukov, Cinelli, Newey, Sharma and Syrgkanis (2026) covers it, "
"and the bound is well posed here."
)
#: Why the plug-in estimator of nu^2 publishes no confidence limit. One text for the six
#: guarded accessors of :class:`SensitivityBounds`, its summary line and the robustness row
#: of the assessment, which all read it through :func:`limit_refusal_reason`.
#: ``tests/unit/test_documentation_links.py`` checks that it cites F26.
_PLUGIN_LIMITS_REFUSAL = (
"no derivation in a source this package cites gives the standard error of a bound built "
"on the plug-in nu^2. E_n[alpha_hat^2] squares the fitted Riesz representer, so it moves "
"at first order with the fitted treatment mechanism, and its curve alpha_hat^2 - nu^2 has "
"no term for that fit. Lemma 3 and Theorem 4 of Chernozhukov, Cinelli, Newey, Sharma and "
"Syrgkanis (2026) give the limits for the doubly robust estimator only. The bias-adjusted "
"bounds, rv, max_bias, benchmark() and contour() remain available, and "
"nu2_estimator='doubly_robust' gives the limits for any mechanism learner at the rates "
"of that Theorem 4. "
"docs/roadmap.md F26 tracks this stop, and docs/roadmap.md X26 tracks the limits for a "
"fit whose own treatment learner is a declared unpenalized logistic or multinomial model"
)
#: The estimators of nu^2 whose bounds refuse their limits, and why. The doubly robust
#: estimator is absent: its limits are the ones Theorem 4 derives.
_ESTIMATOR_LIMIT_REFUSALS: dict[str, str] = {
"plugin": _PLUGIN_LIMITS_REFUSAL,
}
_ZERO_BIAS_LIMITS_REFUSAL = (
"max_bias is zero: sqrt(sigma^2 nu^2) has no ordinary influence curve at "
"sigma^2 = 0, so the one-sided limits and confidence-limit robustness value "
"are not available"
)
_UNREACHABLE_RV = "no equal strength reaches the null"
def limit_refusal_reason(nu2_estimator: str) -> str | None:
"""Why a bound built by ``nu2_estimator`` publishes no limit, or ``None`` when it does.
The accessors of :class:`SensitivityBounds`, its summary and the robustness row of the
assessment each prefix their own subject to this one text, so the three cannot disagree.
Parameters
----------
nu2_estimator : str
The resolved estimator a bound or a robustness report records.
Returns
-------
str or None
``"nu2_estimator='plugin': <reason>."``, or ``None`` for the doubly robust estimator.
"""
reason = _ESTIMATOR_LIMIT_REFUSALS.get(nu2_estimator)
if reason is None:
return None
return f"nu2_estimator={nu2_estimator!r}: {reason}."
_RESPONSE_BOUND_REFUSAL = (
"the omitted-variable bound is not implemented for a fit with a response mechanism. "
"The identified mean is a linear functional of the regression of Delta Y on "
"(A, Delta, W), " + _THEOREM_2_COVERS + " Three pieces are missing. The implemented "
"representer omits the response indicator Delta. sigma^2 averages the respondents, "
"and the theorem needs E[Delta (Y - Qbar)^2] over every row. cf_d would measure a "
"joint strength over the treatment and response mechanisms, and this package reports "
"it as a treatment-side strength."
)
_INTERMEDIATE_BOUND_REFUSAL = Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"the omitted-variable bound is not implemented for a fit with an intermediate "
"variable. Each estimand at the fixed level z is a linear functional of the "
"regression of Y on (A, Z, W), " + _THEOREM_2_COVERS + " The representer carries the "
"intermediate weight 1{Z = z} / P(Z = z | A, W), so cf_d would measure a joint "
"strength over the treatment and intermediate mechanisms. This package reports cf_d "
"and the robustness value as a treatment-side strength.",
alternative=NO_ALTERNATIVE,
roadmap="X9",
)
#: Why the median of several fold draws has no bound, in the words
#: :func:`~cleverly.exceptions.repeats_refusal` completes.
_REPEATS_REASON = (
"A coordinatewise median of the bound's influence terms would not be the influence "
"function of the median bound."
)
#: Appended to the response refusal on a fit that can still run the tilt. A
#: natural-course, shift, incremental, regime, MSM, or ratio-only fit cannot: both tilt
#: rows are unavailable there
#: (:func:`~cleverly.sensitivity.missingness.fit_wide_tilt_refusal`), so the pointer would
#: send the reader to a second refusal.
_RESPONSE_TILT_POINTER = (
" The missingness tilt remains the sensitivity analysis for response: call "
"sensitivity.missingness() or sensitivity.tipping_gamma()."
)
#: The pointer on a fit that also reports a target read from the natural-course mean.
#: The tilt covers the arm means only, so the pointer says so rather than sending a PAR
#: or PAF request to a tilt that refuses it.
_RESPONSE_ARM_TILT_POINTER = (
" The missingness tilt remains the sensitivity analysis for response on this fit's "
"arm means, and not on a target read from the natural-course mean: call "
"sensitivity.missingness() or sensitivity.tipping_gamma() and name the arm means."
)
def _refuse_longitudinal(result: Any) -> Refusal | None:
"""Refuse any assessment family but ``point``.
First in the table, and not by taste: :class:`~cleverly.longitudinal.LongitudinalData`
declares no ``has_missing_outcome``, so a later rule would raise ``AttributeError`` on
a longitudinal result instead of refusing it.
"""
if getattr(result, "assessment_family", None) != "point":
return _LONGITUDINAL_BOUND_REFUSAL
return None
def _refuse_guarded_mechanism(result: Any) -> Refusal | None:
"""Refuse a DR-TMLE fit, whose estimator does not assume a consistent mechanism."""
if result.fitted_method == "drtmle":
return _DRTMLE_BOUND_REFUSAL
return None
def _refuse_selected_mechanism(result: Any) -> Refusal | None:
"""Refuse a collaborative fit, whose representer comes from a selected working g.
The fits that select nothing have their own sentence. The selection records the
declared ``"discrete"`` list as its path, and a ``"learners"`` path is the declared list
followed by the terminal candidate, so the path without its last entry stands for the
declared ``"learners"`` list.
A ``"derivative"`` fit (C-TMLE0) selects nothing and assumes a consistent mechanism as
the ordinary TMLE does, so the bound admits it. Its ``benchmark`` refit meets the refit
refusal of a fit that declared propensity arrays.
"""
if result.fitted_method != "collaborative_tmle":
return None
selection = result.extra.get("ctmle")
if isinstance(selection, CTMLEDerivativeFit):
return None
if isinstance(selection, CTMLESelection) and declares_no_selection(
selection.strategy,
selection.path,
selection.path[:-1] if selection.strategy == "learners" else None,
selection.covariates,
):
return _CTMLE_FULL_CANDIDATE_BOUND_REFUSAL
return _CTMLE_BOUND_REFUSAL
def _refuse_response_mechanism(result: Any) -> Refusal | None:
"""Refuse a fit whose outcome is unobserved on some rows.
The predicate is the data flag ``has_missing_outcome``, which the missingness-tilt
table reads too, so every fit those rows offer the tilt on is a fit this rule refuses.
It reads the data rather than a fitted missingness nuisance on one repeat, because a
response mechanism is a property of the identified functional and survives a
replacement of the stored nuisances. The sibling surface
:mod:`~cleverly.sensitivity._simulated_confounding_request` reads the same flag for
the same boundary. ``result.data`` has no default here, so a result without
point-treatment data raises ``AttributeError`` rather than reporting no mechanism.
The sentence points at the tilt only when
:func:`~cleverly.sensitivity.missingness.fit_wide_tilt_refusal` admits the fit, which
is the predicate both tilt rows read. Any narrower test, such as a natural-course
check alone, would send a shift or a regime fit to a tilt that refuses.
"""
if not result.data.has_missing_outcome:
return None
if fit_wide_tilt_refusal(result) is not None:
pointer = ""
elif reads_natural_course_mean(result):
pointer = _RESPONSE_ARM_TILT_POINTER
else:
pointer = _RESPONSE_TILT_POINTER
return Refusal(
RefusalReason.IMPLEMENTATION_GAP,
_RESPONSE_BOUND_REFUSAL,
alternative=pointer.strip() or NO_ALTERNATIVE,
roadmap="X9",
)
#: Why the bound refuses a fit with a declared missing treatment.
_MISSING_TREATMENT_BOUND_REFUSAL = Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"the omitted-variable bound is not implemented for a fit with a declared missing "
"treatment. The Riesz representer of the bound is derived for a recorded treatment, "
"and this fit divides by the composite mechanism P(A = a, Delta_A = 1[, Delta = 1] | W), "
"whose representer and confounding strengths are not derived. Roadmap X9 owns the "
"representers that carry a second mechanism.",
alternative="Call evalue(result) on the same fit, which has no missing-treatment rule.",
roadmap="X9",
)
def _refuse_missing_treatment(result: Any) -> Refusal | None:
"""Refuse a fit with a declared missing treatment, whose representer is not derived."""
if not getattr(result.data, "has_missing_treatment", False):
return None
return _MISSING_TREATMENT_BOUND_REFUSAL
def _refuse_intermediate(result: Any) -> Refusal | None:
"""Refuse a fit that declares an intermediate variable.
The predicate is :func:`~cleverly.estimators.direct_effect.declares_intermediate`,
which every sensitivity surface reads for this boundary, so they refuse the same fits.
Without this rule the bound returns a number, and nothing in it says that ``cf_d`` has
changed meaning.
"""
return _INTERMEDIATE_BOUND_REFUSAL if declares_intermediate(result) else None
def _refuse_non_arm_axis(result: Any) -> Refusal | None:
"""Refuse a fit whose counterfactuals are indexed by anything but an arm.
Two reasons, not one. A regime, a shift, and an MSM coefficient are linear
functionals of the outcome regression and each has a Riesz representer, so the bound
is well posed and unimplemented. An incremental intervention builds its density out
of the estimated mechanism, so the mechanism is part of the estimand rather than a
nuisance the bound conditions on, and the regression-only bound does not cover it.
"""
axis = result.config.parameter_axis
if axis == "arm":
return None
if axis == "learned_rule":
# A third reason: the target is data-adaptive, and no bound for it was reviewed.
# A benchmark refit would also relearn the rules, and so change the target.
return Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"the omitted-variable bound and its benchmark do not cover a learned-rule fit: "
+ LEARNED_RULE_BOUND_REFUSAL
+ ".",
alternative=NO_ALTERNATIVE,
roadmap="X11",
)
if axis == "ipsi":
return Refusal(
RefusalReason.INVALID_DECLARATION,
"the omitted-variable bound does not cover a fit whose parameters are indexed "
"by 'ipsi'. An incremental intervention tilts the treatment mechanism, so the "
"mechanism is part of the estimand rather than a nuisance the bound conditions "
"on. Chernozhukov, Cinelli, Newey, Sharma and Syrgkanis (2026) bound the bias "
"of a linear functional of the outcome regression alone.",
)
noun = _AXIS_NOUNS[axis]
return Refusal(
RefusalReason.IMPLEMENTATION_GAP,
f"the omitted-variable bound is not implemented for a fit whose parameters are "
f"indexed by {axis!r}. {noun} is a linear functional of the outcome regression, "
f"and it has a Riesz representer, so the bound is well posed here and no "
f"implementation of it is registered. The implemented bound covers the "
f"arm-indexed linear estimands {sorted(LINEAR_ESTIMANDS)}.",
alternative=NO_ALTERNATIVE,
roadmap="X9",
)
def _refuse_unbounded_parameters(result: Any) -> Refusal | None:
"""Refuse an arm-indexed fit that reports no parameter the bound applies to.
The predicate is :func:`~cleverly.sensitivity._parameters.reported_arm_parameters`,
which :func:`resolve_parameter` reads for the same fit, so no request can lift this
rule. A ratio-only fit, a population attributable risk or fraction, and a
complete-outcome natural-course mean each report none, so every call on such a fit
refuses, and the five rows read this rule. The rule follows ``parameter_axis``,
because its sentence describes an arm-indexed fit. The E-value pointer is appended
only on a fit that reports a risk ratio or an odds ratio, the two scales the E-value
answers for.
"""
if reported_arm_parameters(result):
return None
ratio = any(parameter_stem(name) in _EVALUE_STEMS for name in result.estimates)
return Refusal(
RefusalReason.INVALID_DECLARATION,
"the omitted-variable bound applies to the arm-indexed linear estimands "
f"{sorted(LINEAR_ESTIMANDS)}, and this arm-indexed fit reported none of them: "
f"it reported {sorted(result.estimates)}. One bound is the second moment of "
"one contrast's own Riesz representer, so ask the fit for a counterfactual "
"mean or a contrast of arms.",
alternative=_EVALUE_POINTER.strip() if ratio else None,
)
def _refuse_repeats(result: Any) -> Refusal | None:
"""Refuse a median-combined repeated fit.
The one rule a refit of the same estimator lifts, which is why it runs last. Every
rule above it refuses the fit at any split count, so a reader told to fit one split
meets no second refusal after the refit.
"""
return repeats_refusal(
result.n_repeats,
operation="omitted-variable sensitivity",
text=_REPEATS_REASON,
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative="Fit one split for this analysis.",
roadmap="X65",
)
#: Every omitted-variable refusal that the requested estimand cannot change, in the one
#: order the five entry points and the capability rows both use. The table is ordered and
#: each rule assumes its predecessors returned ``None``: ``longitudinal`` establishes that
#: the result carries a point-treatment ``data`` and ``config`` at all, which every rule
#: after it reads. ``bound_parameters`` follows ``parameter_axis``, because its sentence
#: names an arm-indexed fit. ``repeats`` is last because it is the only rule a refit with
#: one split lifts. The names are the introspection contract; a test reads them to pin the
#: order without respelling a message.
_FIT_WIDE_BOUND_RULES: tuple[tuple[str, Callable[[Any], Refusal | None]], ...] = (
("longitudinal", _refuse_longitudinal),
("drtmle", _refuse_guarded_mechanism),
("collaborative_tmle", _refuse_selected_mechanism),
("response_mechanism", _refuse_response_mechanism),
("missing_treatment", _refuse_missing_treatment),
("intermediate", _refuse_intermediate),
("parameter_axis", _refuse_non_arm_axis),
("bound_parameters", _refuse_unbounded_parameters),
("repeats", _refuse_repeats),
)
def fit_wide_bound_refusal(result: Any) -> Refusal | None:
"""Return the first omitted-variable refusal that applies to the whole fit.
Every boundary in :data:`_FIT_WIDE_BOUND_RULES` is reachable from capability reporting
and from execution, and the requested estimand cannot change its verdict. One helper
answers both callers, so a fit the bound refuses is never advertised as available.
Estimand-specific checks stay in :func:`resolve_parameter`.
Parameters
----------
result : Any
Fitted result inspected by this module or by its assessment facade.
Returns
-------
Refusal or None
The first fit-wide refusal, or ``None`` when no fit-wide boundary applies.
"""
for _name, rule in _FIT_WIDE_BOUND_RULES:
refusal = rule(result)
if refusal is not None:
return refusal
return None
[docs]
@dataclass(frozen=True)
class SensitivityElements:
r"""The ingredients of the bias bound for one estimand.
Parameters
----------
estimand : str
Alias of the estimand these elements describe.
sigma2 : float
:math:`E[(Y - \bar Q^*(A, W))^2]`, the residual outcome variance.
nu2 : float
:math:`E[\alpha(A, W)^2]`, the second moment of the Riesz representer.
max_bias : float
:math:`\sqrt{\sigma^2 \nu^2}`, the bias scale. The bound multiplies it by
``|rho| sqrt(cf_y cf_d / (1 - cf_d))``, which can exceed 1, so it is not the
largest bias. The name follows DoubleML.
psi_sigma2 : ndarray
Influence curve of ``sigma2``.
psi_nu2 : ndarray or None
Influence curve of ``nu2`` under the doubly robust estimator: the score of Lemma 3
of Chernozhukov et al. (2026), plus the influence of the conditioning share for the
ATT and the ATC. ``None`` under the plug-in estimator, whose curve no derivation
in a source this package cites gives (docs/roadmap.md F26).
psi_max_bias : ndarray or None
Influence curve of ``max_bias``, as in Theorem 4 of the same paper, so the
bias-adjusted bounds get confidence limits rather than being treated as known
constants. ``None`` under the plug-in estimator or when ``max_bias`` is zero.
riesz_representer : ndarray
``(n,)`` values of :math:`\\alpha(A, W)` for the targeted functional.
nu2_estimator : str
Which estimator of ``nu2`` produced these values.
"""
estimand: str
sigma2: float
nu2: float
max_bias: float
psi_sigma2: FloatArray
psi_nu2: FloatArray | None
psi_max_bias: FloatArray | None
riesz_representer: FloatArray
nu2_estimator: str
def sensitivity_elements(
result: TMLEResult,
estimand: str = "ate",
*,
nu2_estimator: str = "auto",
) -> SensitivityElements:
r"""Compute :math:`\sigma^2`, :math:`\nu^2` and the bias scale for one estimand.
This is where every entry point in this module meets
:data:`_FIT_WIDE_BOUND_RULES`, so :func:`omitted_variable_bounds`,
:func:`benchmark`, :func:`robustness_value`, and :func:`contour_data` share one
refusal and none of them pays a learner before hearing it.
Parameters
----------
result : TMLEResult
A fitted result.
estimand : str
Alias to build the elements for.
nu2_estimator : {"auto", "doubly_robust", "plugin"}
``"doubly_robust"`` uses :math:`E[2 m(O, \alpha) - \alpha^2]`, which is less
sensitive to error in the estimated propensity than the plug-in
:math:`E[\alpha^2]` (``"plugin"``). Both are consistent; ``"auto"`` picks the
doubly robust form wherever the functional's :math:`m(O, \alpha)` has a closed
form, which is all of :data:`LINEAR_ESTIMANDS`. Only the doubly robust form has
a derived influence curve, so only it gives the bounds confidence limits.
Returns
-------
SensitivityElements
The residual outcome variance, the Riesz second moment, the bias scale, and
the influence curve of each.
Raises
------
ValueError
If ``nu2_estimator`` is not one of :data:`NU2_ESTIMATORS`. The check runs
before any refusal and before any computation.
CapabilityError
If a rule of :data:`_FIT_WIDE_BOUND_RULES` refuses the fit, if the fit reports
no parameter the bound applies to under ``estimand``, or if the doubly robust
estimate of :math:`\nu^2` is not positive.
"""
_resolved_nu2_estimator(nu2_estimator)
refusal = fit_wide_bound_refusal(result)
if refusal is not None:
# The tilt pointer is fit-wide: a joint natural-course fit can tilt its arm
# means. The tilt refuses a target read from the natural-course mean, so a
# request for one is not sent there.
if natural_course_names([estimand]) and refusal.alternative in {
_RESPONSE_TILT_POINTER.strip(),
_RESPONSE_ARM_TILT_POINTER.strip(),
}:
refusal = replace(refusal, alternative=NO_ALTERNATIVE)
raise CapabilityError.from_refusal(refusal)
parameter = resolve_parameter(result, estimand)
return _elements_for(result, result.repeats[0], parameter, nu2_estimator)
def _resolved_nu2_estimator(nu2_estimator: str) -> str:
"""The estimator ``nu2_estimator`` names, with ``"auto"`` resolved, or a ``ValueError``."""
if nu2_estimator not in NU2_ESTIMATORS:
raise ValueError(
f"nu2_estimator must be one of {list(NU2_ESTIMATORS)}; got {nu2_estimator!r}"
)
return "doubly_robust" if nu2_estimator == "auto" else nu2_estimator
def resolve_parameter(result: TMLEResult, estimand: str) -> ArmParameter:
"""The arms a requested estimand is about, or a refusal that says why not.
One bound is one linear functional, so this is where "which contrast" is decided --
``ate`` on a two-armed fit and ``ate[medium vs low]`` on a wider one, each with its
own Riesz representer. It reads the request alone. The caller has already passed
the fit through :func:`fit_wide_bound_refusal`, whose ``parameter_axis`` rule makes
the fit arm-indexed and whose ``bound_parameters`` rule makes it report at least one
parameter the bound applies to, since no estimand name can change either verdict.
The order of the checks below matters. A name this bound could never apply to is
refused first, so that asking for a risk ratio is not reported as a missing
estimand; the coverage message comes last.
Parameters
----------
result : TMLEResult
A fitted result.
estimand : str
Alias the caller asked for.
Returns
-------
ArmParameter
The reported parameter and the arms it is a functional of.
"""
stem = parameter_stem(estimand)
if stem not in LINEAR_ESTIMANDS:
# The pointer is for the two ratio scales an E-value is defined on. Offered for
# anything else -- an attributable fraction, say -- it sends the reader from one
# refusal to a second one.
pointer = _EVALUE_POINTER.strip() if stem in _EVALUE_STEMS else None
raise CapabilityError(
f"the omitted-variable bound applies to {sorted(LINEAR_ESTIMANDS)}, not "
f"{estimand!r}: the bound is on the bias of a linear functional of the "
f"outcome regression, and {estimand!r} is not one.",
reason=RefusalReason.INVALID_DECLARATION,
alternative=pointer,
)
available = reported_arm_parameters(result)
if estimand in available:
return available[estimand]
# Before the coverage message: this one *was* reported, so "not requested" would be
# false. It is the derivation that is missing, not the parameter.
conditional = stratum_refusal(result, estimand, "the omitted-variable bound")
if conditional is not None:
raise CapabilityError.from_refusal(conditional)
empty = _refuse_unbounded_parameters(result)
if empty is not None:
# Every entry point hears this sentence from the fit-wide table before it gets
# here. A caller that resolves a parameter directly hears the same sentence.
raise CapabilityError.from_refusal(empty)
raise CapabilityError(
f"estimand {estimand!r} was not requested in this fit. The bound is available for "
f"{sorted(available)}, one per contrast, since nu^2 is the second moment of "
"that contrast's own Riesz representer.",
reason=RefusalReason.INVALID_DECLARATION,
)
def _elements_for(
result: TMLEResult,
repeat: RepeatFit,
parameter: ArmParameter,
nu2_estimator: str,
) -> SensitivityElements:
"""The bound's pieces under one cross-fitting draw.
Takes the targeted ``Qbar`` and the mechanism from the same
:class:`~cleverly.estimators._nuisance.RepeatFit`, which is what makes ``sigma2`` the
residual variance of the regression whose propensity ``nu2`` was computed from.
"""
data = result.data
scaler = repeat.nuisance.scaler
group = parameter.group
fluctuation = repeat.fluctuations[group]
bounds = g_bounds_for(group, result.config.g_bounds, result.config.g_bounds_conditional)
reference = result.config.reference_arm
submodel = build_submodel(
data,
repeat.nuisance,
group,
bounds=bounds,
nuisance_bound=result.config.missingness_bound,
intermediate_value=result.intermediate_value,
# The conditional-effect fluctuations contrast against the arm this fit declared,
# so the covariate rebuilt here must be the one it was targeted with.
reference=reference,
)
# The weight of the arm the estimand *conditions on*: ``_m_alpha`` multiplies the
# contrast by 1{A = c} / P(A = c), the observed membership of that arm over its share,
# and never by a fitted g_c. ``None`` for a mean or an unconditional contrast, which
# reweight nobody -- and the arm is read off the parameter rather than assumed to be
# the other one, since with K arms ``att[medium vs low]`` and ``att[high vs low]``
# condition on different populations.
conditioning = parameter.conditions_on
conditioning_indicator: FloatArray | None = None
conditioning_share: float | None = None
if conditioning is not None:
conditioning_indicator = np.asarray(data.treatment == conditioning, dtype=float)
conditioning_share = float(data.arm_fractions[repeat.nuisance.arms.index(conditioning)])
# sigma^2: residual variance of the targeted outcome regression, on the original
# outcome scale so the bound is reported in the units the estimate uses.
scaled = scaler.scale(data.outcome)
residual = np.where(data.observed, scaled - fluctuation.targeted.observed, 0.0)
if not scaler.is_identity:
residual = residual * scaler.range
weights = data.weights
sigma2_element = residual**2
sigma2 = float(np.average(sigma2_element[data.observed], weights=weights[data.observed]))
psi_sigma2 = np.where(data.observed, sigma2_element - sigma2, 0.0) * weights
representer = _riesz_representer(parameter, submodel)
method = _resolved_nu2_estimator(nu2_estimator)
if method == "doubly_robust":
m_alpha = _m_alpha(parameter, submodel, conditioning_indicator, conditioning_share)
nu2_element = 2.0 * m_alpha - representer**2
nu2 = float(np.average(nu2_element, weights=weights))
if nu2 <= 0:
# A refusal rather than the plug-in value. By the Riesz identity this
# estimator is nu_0^2 minus the squared error of the fitted representer, so a
# nonpositive value reports a mechanism fit further from the truth than the
# representer is large. The plug-in E[alpha_hat^2] is built from that same
# fitted representer, so substituting it reports a second moment of the wrong
# function under the name of the right one.
# ``CapabilityError`` rather than ``ValueError``: ``run_all`` catches only the
# refusal type, and a bare ``ValueError`` would abort the whole battery.
raise CapabilityError(
f"the {method!r} estimator of nu^2 returned {nu2:.6g} for "
f"{parameter.name!r}, and a second moment cannot be negative. "
f"E[2 m(alpha_hat) - alpha_hat^2] equals nu_0^2 minus the squared error of "
f"the fitted representer, so this value reports a treatment mechanism the "
f"bound's derivation does not cover. The plug-in E[alpha_hat^2] squares "
f"that same fitted representer, so it is not a substitute. "
f"(nu2_estimator={nu2_estimator!r} resolved to {method!r}.)",
reason=RefusalReason.INVALID_DECLARATION,
)
else:
nu2_element = representer**2
nu2 = float(np.average(nu2_element, weights=weights))
max_bias = float(np.sqrt(sigma2 * nu2))
psi_nu2: FloatArray | None = None
psi_max_bias: FloatArray | None = None
# The plug-in keeps no curve: E_n[alpha_hat^2] moves at first order with the fitted
# mechanism and no derivation in a source this package cites gives its influence
# function (docs/roadmap.md F26), so an
# array here would be a curve for a different estimator under this one's name.
if method == "doubly_robust":
psi_nu2 = (nu2_element - nu2) * weights
if conditioning_indicator is not None:
# The estimate divides by the estimated share of the conditioning arm, so its
# curve carries that share's influence. See _conditioning_share_influence.
assert conditioning_share is not None
psi_nu2 = psi_nu2 + _conditioning_share_influence(
nu2, conditioning_indicator, conditioning_share, weights
)
# The square-root delta method requires a positive product. At sigma2=0 the
# point bound is still defined, but its ordinary influence curve is not.
if max_bias > 0.0:
psi_max_bias = (sigma2 * psi_nu2 + nu2 * psi_sigma2) / (2.0 * max_bias)
return SensitivityElements(
estimand=parameter.name,
sigma2=sigma2,
nu2=nu2,
max_bias=max_bias,
psi_sigma2=psi_sigma2,
psi_nu2=psi_nu2,
psi_max_bias=psi_max_bias,
riesz_representer=representer,
nu2_estimator=method,
)
def _riesz_representer(parameter: ArmParameter, submodel: Any) -> FloatArray:
r"""``alpha(A, W)`` for the requested estimand.
The clever covariate *is* the Riesz representer -- that is why the same object
both drives the targeting step and controls how much an omitted confounder can
move the estimate. Poor overlap inflates both simultaneously.
Every column is reached by the arm it carries rather than by position, so
``ate[high vs low]`` reads the two arms it names out of a ``K``-column covariate
instead of the two an implementation that counted to two would have found.
"""
if parameter.group in ("att", "atc"):
# One column per non-reference arm, and this parameter's is the one carrying the
# contrast it is named for.
return submodel.contrast_column_for(parameter.arm)
if parameter.versus is None:
return submodel.column_for(parameter.arm)
return np.asarray(
submodel.column_for(parameter.arm) - submodel.column_for(parameter.versus), dtype=float
)
def _m_alpha(
parameter: ArmParameter,
submodel: Any,
conditioning_indicator: FloatArray | None,
conditioning_share: float | None,
) -> FloatArray:
r"""The target functional applied to the Riesz representer, ``m(O, alpha)``.
Used by the doubly robust estimator of :math:`\nu^2`, which relies on the Riesz
identity :math:`E[m(O, \alpha)] = E[\alpha_0 \alpha]` for every :math:`\alpha`. Put
:math:`\hat\alpha` in it and :math:`E[2 m(O, \hat\alpha) - \hat\alpha^2]` equals
:math:`\nu_0^2 - E[(\hat\alpha - \alpha_0)^2]`, which is what makes the estimator
first-order insensitive to the fitted mechanism.
The identity holds only when ``m`` is the functional's own score, with no fitted
nuisance in it. An ATT averages the contrast over the units that received the
conditioning arm ``c``, so its score weights the contrast by the observed
:math:`1\{A = c\} / P(A = c)`. That is the ATT of Example 2 in Online Appendix A of
Chernozhukov, Cinelli, Newey, Sharma and Syrgkanis (2026), with
:math:`\ell(W_s) = D / P(D = 1)`. Item (2) of their Theorem 5 writes this :math:`m`
and the representer. That theorem establishes the bias formula; the score and the
limits are Lemma 3 and Theorem 4. DoubleML's
``DoubleMLIRM._sensitivity_element_est`` writes the same score, expanded, as
:math:`D / (p^2 (1 - \hat m))`. A fitted :math:`\hat g_c(W) / P(A = c)` in that
place has the right mean only at :math:`\hat g = g_0`. Anywhere else the identity
fails, the estimator loses its insensitivity to the mechanism, and it can exceed
:math:`\nu_0^2`.
``conditioning_indicator`` and ``conditioning_share`` belong to the arm the estimand
conditions on -- the contrast arm for an ATT and the reference for an ATC -- rather
than to arm 1 and its complement. With two arms and the default reference those are
the same two numbers; naming the arm is what keeps them right when a binary fit
declares the other one, and what gives each of ``K - 1`` contrasts its own
population. They are ``None`` for a mean or an unconditional contrast, which
reweight nobody.
"""
# ``arms[a][:, c]`` is the covariate at arm ``a`` in the column targeting arm ``c``:
# the mean submodel has one column per arm, so both indices are arm levels and
# neither is a positional 0 or 1.
if parameter.group == "mean":
columns = submodel.arm_columns
at = parameter.arm
if parameter.versus is None:
return np.asarray(submodel.arms[at][:, columns[at]], dtype=float)
versus = parameter.versus
first, second = submodel.arms[at], submodel.arms[versus]
return np.asarray(
first[:, columns[at]]
- first[:, columns[versus]]
- (second[:, columns[at]] - second[:, columns[versus]]),
dtype=float,
)
# ATT / ATC: the functional carries the observed arm-membership weight, since it
# averages the contrast over the conditioning arm's subpopulation rather than over
# everyone. The column is read by the arm whose contrast it carries, not as a literal 0.
assert parameter.versus is not None
assert conditioning_indicator is not None and conditioning_share is not None
column = submodel.contrast_columns[parameter.arm]
difference = np.asarray(
submodel.arms[parameter.arm][:, column] - submodel.arms[parameter.versus][:, column],
dtype=float,
)
return np.asarray((conditioning_indicator / conditioning_share) * difference, dtype=float)
def _conditioning_share_influence(
nu2: float, indicator: FloatArray, share: float, weights: FloatArray
) -> FloatArray:
r"""The influence of the estimated conditioning share on the doubly robust ``nu^2``.
For the ATT and the ATC, the representer and :math:`m(O, \alpha)` both divide by the
share :math:`p = P_w(A = c)` of the arm :math:`c` that the parameter conditions on. So
the estimate is :math:`\hat F / \hat p^2`, where :math:`\hat F` is the weighted mean of
:math:`\hat p^2 (2 m(O, \hat\alpha) - \hat\alpha^2)`. At a fixed :math:`p`, the Riesz
identity makes :math:`\hat F` insensitive at first order to the fitted mechanism, so its
curve is the score of Lemma 3 of Chernozhukov, Cinelli, Newey, Sharma and Syrgkanis
(2026). The weighted share :math:`\hat p` has the curve :math:`w (1\{A = c\} - p)`, and
the derivative of :math:`F / p^2` in :math:`p` is :math:`-2 \nu^2 / p`. The delta method
therefore adds
.. math::
-2 \nu^2 w (1\{A = c\} - p) / p
to the curve. With ``K`` arms each contrast has its own :math:`c` and :math:`p`. Under
clustering the term enters the cluster sum of
:func:`~cleverly.inference.cluster.influence_variance` like every other row term.
The paper gives a term of this kind for the estimate only, in Online Appendix A,
"Statistical Inference", Equation (15), and the estimate's own curve here carries it.
It gives none for :math:`\nu^2`, so this is the package's derivation. The R package
``dml.sensemakr`` carries the same term since commit ``d5293ecb``, and DoubleML's
``DoubleMLIRM._sensitivity_element_est`` omits it. The weighted rows of the term sum
to zero, so no check of a mean or of a point can see it.
``tests/unit/test_omitted_variable_standard_error.py`` compares the curve with its
Gateaux derivative row by row.
"""
return np.asarray(-2.0 * nu2 * (indicator - share) / share * weights, dtype=float)
_P = ParamSpec("_P")
#: The constructor arguments whose values the bound stores under a private field name.
_STORED_LIMITS: tuple[str, ...] = ("ci_lower", "ci_upper", "robustness_value_ci")
def _accepts_stored_limit_names(init: Callable[_P, None]) -> Callable[_P, None]:
"""Let the constructor take the private field names, so that ``replace()`` works.
``dataclasses.replace()`` and ``copy.replace()`` pass every field under its stored
name, and :class:`SensitivityBounds` stores its three limits as ``_ci_lower``,
``_ci_upper`` and ``_robustness_value_ci`` so that the public names can refuse. The
wrapper renames those keywords to the public arguments, and a public argument passed
beside its stored name wins, as in ``replace(bound, ci_lower=x)``. ``functools.wraps``
keeps the public signature for :func:`inspect.signature`, and the wrapped constructor
raises its own ``TypeError`` on an unknown or a missing argument.
"""
@functools.wraps(init)
def wrapper(*args: _P.args, **kwargs: _P.kwargs) -> None:
for name in _STORED_LIMITS:
if f"_{name}" in kwargs:
stored = kwargs.pop(f"_{name}")
kwargs.setdefault(name, stored)
init(*args, **kwargs)
return wrapper
[docs]
@dataclass(frozen=True, init=False)
class SensitivityBounds:
"""Bias-adjusted bounds under an assumed confounder strength.
Parameters
----------
estimand : str
Alias of the estimand these bounds describe.
psi : float
The unadjusted point estimate.
cf_y : float
Share of the residual outcome variation the assumed confounder explains.
cf_d : float
Share of the Riesz representer's second moment that only the assumed confounder
adds. For the ATE it is the gain in precision of the treatment model.
rho : float
How adversarially the two are aligned. ``1.0`` is the worst case.
confounding_strength : float
The product those three imply, ``|rho| sqrt(cf_y cf_d / (1 - cf_d))``.
max_bias : float
The bias scale ``sqrt(sigma^2 nu^2)``. :attr:`bias` multiplies it by
``confounding_strength``, which can exceed 1. The name follows DoubleML.
lower : float
Bias-adjusted lower bound on the estimate.
upper : float
Bias-adjusted upper bound on the estimate.
ci_lower : float
Lower one-sided confidence limit of the adjusted bound. It refuses at a
non-inferential status, where :attr:`plugin_interval_lower` keeps the diagnostic.
It refuses under the plug-in estimator of ``nu^2``, which has no derived standard
error (docs/roadmap.md F26).
ci_upper : float
Upper one-sided confidence limit of the adjusted bound. It follows ``ci_lower``.
level : float
Coverage level of those limits.
robustness_value : float or None
Equal confounding strength that moves the point estimate to the null. ``None``
means no strength below one reaches it.
robustness_value_ci : float or None
The same strength for the confidence limit. It refuses with ``ci_lower``;
:attr:`robustness_value_plugin_interval` keeps the diagnostic. ``None`` means
the available limit does not reach the null at any equal strength below one.
null_hypothesis : float
The value the robustness values are measured against.
inference : str, default="influence_curve"
The inference status of the estimate the bound adjusts. One of
:data:`~cleverly.inference.influence.InferenceStatus`, which the
:doc:`inference reference </technical-reference/inference>` lists.
nu2_estimator : str, default="doubly_robust"
The resolved estimator of ``nu^2`` that built the bound. Under ``"plugin"`` the six
limit accessors refuse.
"""
estimand: str
psi: float
cf_y: float
cf_d: float
rho: float
confounding_strength: float
max_bias: float
lower: float
upper: float
#: ``None`` when no derivation in a source this package cites gives the limits, so two
#: identical plug-in bounds, and a pickle round trip of one, compare equal. A NaN here
#: compared unequal to itself.
_ci_lower: float | None
_ci_upper: float | None
level: float
robustness_value: float | None
_robustness_value_ci: float | None
null_hypothesis: float
inference: InferenceStatus = "influence_curve"
nu2_estimator: str = "doubly_robust"
@_accepts_stored_limit_names
def __init__(
self,
estimand: str,
psi: float,
cf_y: float,
cf_d: float,
rho: float,
confounding_strength: float,
max_bias: float,
lower: float,
upper: float,
ci_lower: float | None,
ci_upper: float | None,
level: float,
robustness_value: float | None,
robustness_value_ci: float | None,
null_hypothesis: float,
inference: InferenceStatus = "influence_curve",
nu2_estimator: str = "doubly_robust",
) -> None:
# Preserve the public constructor while guarding its inferential values.
for name, value in (
("estimand", estimand),
("psi", psi),
("cf_y", cf_y),
("cf_d", cf_d),
("rho", rho),
("confounding_strength", confounding_strength),
("max_bias", max_bias),
("lower", lower),
("upper", upper),
("_ci_lower", ci_lower),
("_ci_upper", ci_upper),
("level", level),
("robustness_value", robustness_value),
("_robustness_value_ci", robustness_value_ci),
("null_hypothesis", null_hypothesis),
("inference", inference),
("nu2_estimator", nu2_estimator),
):
object.__setattr__(self, name, value)
@property
def _limits_derived(self) -> bool:
"""Whether a derivation gives the standard error of these limits."""
return self.nu2_estimator not in _ESTIMATOR_LIMIT_REFUSALS and self.max_bias > 0.0
def _limit_reason(self) -> str | None:
"""The estimator or zero-bias reason for withholding limit inference."""
reason = limit_refusal_reason(self.nu2_estimator)
if reason is not None:
return reason
if self.max_bias == 0.0:
return _ZERO_BIAS_LIMITS_REFUSAL + "."
return None
def _limit_refusal(self, operation: str, *, diagnostic: bool = False) -> None:
"""Refuse a limit accessor, in the one order every accessor uses.
The inference status comes first, because it is a property of the estimate the
bound adjusts. A diagnostic accessor skips it, since it exists to report the
spread at a non-inferential status. The estimator of ``nu^2`` comes second, and no
name escapes it, because a limit without a derived standard error is not a
diagnostic of anything.
"""
if not diagnostic:
refuse_inference(self.inference, operation=operation)
reason = self._limit_reason()
if reason is not None:
if self.nu2_estimator in _ESTIMATOR_LIMIT_REFUSALS:
record_refusal_site(
"omitted_variable.limit_estimator", RefusalReason.INFERENCE_GAP, "F26"
)
raise CapabilityError(
f"{operation} is not defined under {reason}",
reason=RefusalReason.INFERENCE_GAP,
alternative=(
"Fit the bound with nu2_estimator='doubly_robust', which has "
"derived limits."
),
roadmap="F26",
)
record_refusal_site("omitted_variable.limit_reason", RefusalReason.INFERENCE_GAP, None)
raise CapabilityError(
f"{operation} is not defined under {reason}",
reason=RefusalReason.INFERENCE_GAP,
alternative="Read the point bounds and the robustness value, which stay available.",
)
def _guarded(self, name: str, stored: float | None, *, diagnostic: bool = False) -> float:
"""One stored limit, after :meth:`_limit_refusal` has let it through."""
operation = f"SensitivityBounds.{name}"
self._limit_refusal(operation, diagnostic=diagnostic)
if stored is None: # pragma: no cover - a derived bound always stores its limits
raise CapabilityError(
f"{operation} was not stored on this bound",
reason=RefusalReason.INVALID_DECLARATION,
)
return stored
@property
def ci_lower(self) -> float:
"""Lower confidence limit, when this bound supplies inference."""
return self._guarded("ci_lower", self._ci_lower)
@property
def ci_upper(self) -> float:
"""Upper confidence limit, when this bound supplies inference."""
return self._guarded("ci_upper", self._ci_upper)
@property
def robustness_value_ci(self) -> float | None:
"""Robustness value for a confidence limit, when inference is supplied."""
self._limit_refusal("SensitivityBounds.robustness_value_ci")
return self._robustness_value_ci
@property
def plugin_interval_lower(self) -> float:
"""Lower plug-in limit, a diagnostic at a non-inferential status."""
return self._guarded("plugin_interval_lower", self._ci_lower, diagnostic=True)
@property
def plugin_interval_upper(self) -> float:
"""Upper plug-in limit, a diagnostic at a non-inferential status."""
return self._guarded("plugin_interval_upper", self._ci_upper, diagnostic=True)
@property
def robustness_value_plugin_interval(self) -> float | None:
"""Robustness value from the plug-in limit, a diagnostic when inference is absent."""
self._limit_refusal("SensitivityBounds.robustness_value_plugin_interval", diagnostic=True)
return self._robustness_value_ci
@property
def bias(self) -> float:
"""Magnitude of the bias the assumed confounder could produce."""
return self.confounding_strength * self.max_bias
[docs]
def to_dict(self) -> dict[str, Any]:
"""Return a JSON-compatible representation.
At a non-inferential :attr:`inference` the mapping carries ``inference``, and
the three quantities read off the influence curve take the names
:func:`~cleverly.inference.influence.spread_name` gives them:
``plugin_interval_lower``, ``plugin_interval_upper`` and
``robustness_value_plugin_interval``. When no derivation in a source this package
cites gives the limits' standard error, the mapping carries ``nu2_estimator`` and
omits those three keys under every name.
Returns
-------
dict
A JSON-compatible mapping of every reported field.
"""
status = self.inference
row: dict[str, Any] = {"estimand": self.estimand}
if is_labelled(status):
row["inference"] = status
if not self._limits_derived:
row["nu2_estimator"] = self.nu2_estimator
if self.max_bias == 0.0 and self.nu2_estimator == "doubly_robust":
row["limit_refusal"] = _ZERO_BIAS_LIMITS_REFUSAL
row |= {
"psi": self.psi,
"cf_y": self.cf_y,
"cf_d": self.cf_d,
"rho": self.rho,
"max_bias": self.max_bias,
"bias": self.bias,
"lower": self.lower,
"upper": self.upper,
}
if self._limits_derived:
row |= {
spread_name("ci_lower", status): self.plugin_interval_lower,
spread_name("ci_upper", status): self.plugin_interval_upper,
}
row["robustness_value"] = self.robustness_value
if self._limits_derived:
row[spread_name("robustness_value_ci", status)] = self.robustness_value_plugin_interval
return row
[docs]
def summary(self) -> str:
"""Return a printable summary.
Returns
-------
str
A printable report, one line per reported quantity.
"""
conclusion = (
"the sign of the effect survives"
if (self.lower - self.null_hypothesis) * (self.upper - self.null_hypothesis) > 0
else "the effect could be explained away"
)
header = [
f"Omitted-variable sensitivity for {self.estimand!r}",
"-" * 44,
f"estimate {self.psi:.5g}; bias scale sqrt(sigma^2 nu^2) = {self.max_bias:.5g}",
f"assumed confounding: cf_y = {self.cf_y:.3g}, cf_d = {self.cf_d:.3g}, "
f"rho = {self.rho:.3g}"
f" -> bias <= {self.bias:.5g}",
]
if self.robustness_value is None:
rv_line = f"robustness value RV = unavailable: {_UNREACHABLE_RV}."
else:
rv_line = (
f"robustness value RV = {self.robustness_value:.4f}: a confounder with "
f"cf_y = cf_d = {self.robustness_value:.1%} (the share of the residual outcome "
"variation it explains, and the share of the Riesz representer's second moment "
f"only it adds) would move the estimate to {self.null_hypothesis:g}."
)
reason = self._limit_reason()
if reason is not None:
return "\n".join(
[
*header,
f"bias-adjusted bounds: [{self.lower:.5g}, {self.upper:.5g}] ({conclusion})",
"",
rv_line,
f"no one-sided limits or RVa under {reason}",
]
)
status = self.inference
lower = self.plugin_interval_lower
upper = self.plugin_interval_upper
rv_interval = self.robustness_value_plugin_interval
rva_line = (
f"robustness value RVa = unavailable: {_UNREACHABLE_RV}."
if rv_interval is None
else f"robustness value RVa = {rv_interval:.4f}: the same, for the "
f"{self.level:.0%} {spread_name('confidence bound', status)}."
)
return "\n".join(
[
*header,
f"bias-adjusted bounds: [{self.lower:.5g}, {self.upper:.5g}]",
f"with {self.level:.0%} {spread_name('one-sided CIs', status)}: "
f"[{lower:.5g}, {upper:.5g}] ({conclusion})",
"",
rv_line,
rva_line,
]
)
def __repr__(self) -> str: # pragma: no cover - debugging aid
return self.summary()
def _confounding_strength(cf_y: float, cf_d: float, rho: float) -> float:
"""``|rho| sqrt(cf_y cf_d / (1 - cf_d))``, the multiplier on the bias scale.
``cf_y`` is the share ``eta^2`` of residual outcome variation, so it lies in ``[0, 1]``
and the bound is finite at ``cf_y = 1`` (Theorem 2 of Chernozhukov et al., 2026).
``cf_d = 1`` divides by zero, so ``cf_d`` lies in ``[0, 1)``.
"""
if not 0.0 <= cf_y <= 1.0:
raise ValueError(f"cf_y must lie in [0, 1]; got {cf_y}")
if not 0.0 <= cf_d < 1.0:
raise ValueError(f"cf_d must lie in [0, 1); got {cf_d}. At cf_d = 1 the bound is unbounded")
if not 0.0 <= abs(rho) <= 1.0:
raise ValueError(f"|rho| must lie in [0, 1]; got {rho}")
return float(abs(rho) * np.sqrt(cf_y * cf_d / (1.0 - cf_d)))
def _one_sided_quantile(level: float, df: int | None = None) -> float:
"""Validate the one-sided coverage level before computing its reference quantile.
``df`` is the Student t degrees of freedom of the bounded estimate, or ``None`` for the
normal reference.
"""
if not 0.5 <= level < 1.0:
raise ValueError(f"level must lie in [0.5, 1); got {level}")
return one_sided_quantile(level, df)
def _benchmark_strengths(
result: TMLEResult,
estimand: str,
benchmark: BenchmarkResult,
cf_y: float | None,
cf_d: float | None,
) -> tuple[float, float]:
"""The ``(cf_y, cf_d)`` a benchmark of this fit and estimand assigns to the bound.
The gains enter unchanged. Remark 8 of Appendix E.5 of Chernozhukov et al.
(arXiv:2112.13398v6) states that, with ``k_Y = k_D = 1``, the gain metrics serve as
proxies for the sensitivity parameters themselves. ``share = gain / (1 + gain)`` would
be a different, narrower calibration that no cited source gives.
"""
if not isinstance(benchmark, BenchmarkResult):
raise TypeError(f"benchmark must be a BenchmarkResult; got {type(benchmark).__name__}")
supplied = [name for name, value in (("cf_y", cf_y), ("cf_d", cf_d)) if value is not None]
if supplied:
raise ValueError(
f"benchmark= sets cf_y and cf_d, so {' and '.join(supplied)} cannot also be passed"
)
if benchmark.estimand != estimand:
raise ValueError(
f"the benchmark is of {benchmark.estimand!r}, and this bound is of {estimand!r}. "
f"Pass estimand={benchmark.estimand!r}, or benchmark {estimand!r}"
)
if benchmark.psi_long != result[estimand].psi:
raise ValueError(
f"the benchmark was computed on another fit: its full-fit estimate is "
f"{benchmark.psi_long!r}, and this fit reports {result[estimand].psi!r}"
)
if benchmark.cf_d >= 1.0:
raise ValueError(
f"a confounder as strong as {list(benchmark.covariates)} has a cf_d gain of "
f"{benchmark.gain_d:.4g}. The bound is unbounded at cf_d = 1, so this benchmark "
"calibrates no finite bound"
)
return benchmark.cf_y, benchmark.cf_d
[docs]
def omitted_variable_bounds(
result: TMLEResult,
estimand: str = "ate",
*,
cf_y: float | None = None,
cf_d: float | None = None,
rho: float = 1.0,
level: float = 0.95,
null_hypothesis: float = 0.0,
nu2_estimator: str = "auto",
benchmark: BenchmarkResult | None = None,
) -> SensitivityBounds:
"""Bias-adjusted bounds and robustness values for one estimand.
Defaults follow the convention of assuming a confounder with ``cf_y = cf_d = 0.03``,
with worst-case alignment (``rho = 1``). Prefer
reading :attr:`SensitivityBounds.robustness_value`, which needs no assumption at
all. ``benchmark=`` takes the strengths from a :func:`benchmark` of this fit
instead. Its gain statistics are the published proxies for the two shares
(Chernozhukov et al., arXiv:2112.13398v6, Appendix E.5, Remark 8), so they enter the
bound unchanged.
Parameters
----------
result : TMLEResult
A fitted result.
estimand : str
Alias to bound.
cf_y : float or None
Assumed share of residual outcome variation the confounder explains, in
``[0, 1]``. ``None`` means 0.03, or the benchmark's ``cf_y`` with ``benchmark=``.
cf_d : float or None
Assumed share of the Riesz representer's second moment that only the confounder
adds, in ``[0, 1)``. For the ATE it is the gain in precision of the treatment
model. ``None`` means 0.03, or the benchmark's ``cf_d`` with ``benchmark=``.
rho : float
How adversarially the two are aligned. ``1.0`` is the worst case. A benchmark
does not set it. Pass ``rho=benchmark.rho`` to use the benchmark's alignment.
level : float
Coverage level of the reported limits.
null_hypothesis : float
Value the robustness values are measured against.
nu2_estimator : {"auto", "doubly_robust", "plugin"}
Which estimator of the Riesz second moment to use.
benchmark : BenchmarkResult or None
A benchmark of this fit and this estimand. Its ``cf_y`` and ``cf_d`` become the
assumed strengths, so ``cf_y=`` and ``cf_d=`` must be omitted.
Returns
-------
SensitivityBounds
Adjusted bounds, their confidence limits, and the robustness values. Under
``nu2_estimator="plugin"`` the bounds, ``max_bias`` and the point robustness value
are reported, and the one-sided limits and the confidence-limit robustness value
refuse when read, because no derivation in a source this package cites gives their
standard error (docs/roadmap.md F26). A zero bias scale also leaves the point
bounds available but refuses the limits, whose square-root derivative is undefined.
An unreachable robustness threshold is ``None``.
Raises
------
ValueError
If ``nu2_estimator`` is not one of :data:`NU2_ESTIMATORS`, if ``cf_y`` lies
outside ``[0, 1]`` or ``cf_d`` outside ``[0, 1)``, if ``|rho|`` exceeds one, or if
``level`` lies outside ``[0.5, 1)``. With ``benchmark=``, also if ``cf_y=`` or
``cf_d=`` is passed, or if the benchmark is of another estimand or another fit.
CapabilityError
On every refusal :func:`sensitivity_elements` raises: a fit the rule table
refuses, an estimand the bound does not cover, or a doubly robust estimate of
``nu^2`` that is not positive.
Notes
-----
A fit made with ``treatment_adjustment=`` (``treatment_covariates=`` in
:meth:`cleverly.TMLE.fit`) lets the treatment mechanism read a subset of the adjustment
set. The Riesz representer of the bound uses that fitted mechanism. So the bound reads the
statistical-model restriction that the fit's summary prints, and it is claimed only under
that restriction. Collaborative fits stay refused. If a later release admits a collaborative
fit that selects nothing, the same statement covers it.
"""
_one_sided_quantile(level)
if benchmark is not None:
cf_y, cf_d = _benchmark_strengths(result, estimand, benchmark, cf_y, cf_d)
cf_y = 0.03 if cf_y is None else cf_y
cf_d = 0.03 if cf_d is None else cf_d
elements = sensitivity_elements(result, estimand, nu2_estimator=nu2_estimator)
estimate = result[estimand]
quantile = _one_sided_quantile(level, estimate.reference_df)
strength = _confounding_strength(cf_y, cf_d, rho)
lower = estimate.psi - strength * elements.max_bias
upper = estimate.psi + strength * elements.max_bias
# One-sided limits on each end, Theorem 4 in Section 4 of Chernozhukov et al. (2026).
# The curve adds the uncertainty of the bias term to that of the estimate. The plug-in
# has no curve, so its bound stores no limit and its accessors refuse.
ci_lower: float | None = None
ci_upper: float | None = None
if elements.psi_max_bias is not None:
psi_bias = elements.psi_max_bias
se_lower = _bound_std_error(estimate.influence_curve, -strength * psi_bias, result)
se_upper = _bound_std_error(estimate.influence_curve, strength * psi_bias, result)
ci_lower = lower - quantile * se_lower
ci_upper = upper + quantile * se_upper
rv, rva = _robustness_values(
result, elements, estimate.psi, rho, level, null_hypothesis, estimate.reference_df
)
return SensitivityBounds(
estimand=estimand,
psi=estimate.psi,
cf_y=cf_y,
cf_d=cf_d,
rho=rho,
confounding_strength=strength,
max_bias=elements.max_bias,
lower=lower,
upper=upper,
ci_lower=ci_lower,
ci_upper=ci_upper,
level=level,
robustness_value=rv,
robustness_value_ci=rva,
null_hypothesis=null_hypothesis,
inference=estimate.inference,
nu2_estimator=elements.nu2_estimator,
)
def _bound_std_error(psi_estimate: FloatArray, psi_bias: FloatArray, result: TMLEResult) -> float:
"""Standard error of a bias-adjusted bound."""
combined = np.asarray(psi_estimate, dtype=float) + np.asarray(psi_bias, dtype=float)
return float(np.sqrt(influence_variance(combined, result.data.cluster)))
def _robustness_values(
result: TMLEResult,
elements: SensitivityElements,
psi: float,
rho: float,
level: float,
null_hypothesis: float,
df: int | None = None,
) -> tuple[float | None, float | None]:
"""Solve for ``cf_y = cf_d = v`` at which a bound reaches the null.
The point value is ``None`` when no equal strength in ``[0, 1)`` reaches the null.
The limit value is also ``None`` when its influence curve is unavailable; callers
omit and refuse that accessor rather than interpreting it as an unreachable limit.
"""
side = 1.0 if null_hypothesis > psi else -1.0
psi_bias = elements.psi_max_bias
alignment = abs(rho)
if not 0.0 <= alignment <= 1.0:
raise ValueError(f"|rho| must lie in [0, 1]; got {rho}")
quantile = _one_sided_quantile(level, df)
def equal_share(strength: float) -> float:
assert alignment > 0.0
if strength <= alignment:
ratio = strength / alignment
value = float(2.0 * ratio / (np.hypot(ratio, 2.0) + ratio))
else:
inverse_ratio = alignment / strength
value = float(2.0 / (1.0 + np.sqrt(1.0 + 4.0 * inverse_ratio**2)))
if strength > 0.0 and (not np.isfinite(value) or not 0.0 < value < 1.0):
raise CapabilityError(
"the equal-strength robustness threshold exists but is not representable "
"inside [0, 1) at floating-point precision",
reason=RefusalReason.INVALID_DECLARATION,
)
return value
def limit_at(strength: float) -> float:
assert psi_bias is not None
bias = strength * elements.max_bias
edge = psi + side * bias
se = _bound_std_error(
result[elements.estimand].influence_curve, side * strength * psi_bias, result
)
return edge + side * quantile * se
distance = abs(psi - null_hypothesis)
if distance == 0.0:
rv = 0.0
elif alignment == 0.0 or elements.max_bias == 0.0:
rv = None
else:
rv = equal_share(distance / elements.max_bias)
if psi_bias is None:
return rv, None
# A baseline one-sided limit already on the null side needs no confounding.
# Squared-error minimisation misses this case and reports a small positive value.
def signed_gap(strength: float) -> float:
return side * (limit_at(strength) - null_hypothesis)
if signed_gap(0.0) >= 0.0:
return rv, 0.0
if alignment == 0.0:
return rv, None
high = max(1.0, distance / elements.max_bias)
if not np.isfinite(high):
raise CapabilityError(
"the confidence-limit robustness root exceeds finite precision",
reason=RefusalReason.INVALID_DECLARATION,
)
while signed_gap(high) < 0.0:
if high > np.finfo(float).max / 2.0:
raise CapabilityError(
"the confidence-limit robustness root exceeds finite precision",
reason=RefusalReason.INVALID_DECLARATION,
)
high *= 2.0
strength = float(optimize.brentq(signed_gap, 0.0, high))
return rv, equal_share(strength)
[docs]
@dataclass(frozen=True)
class BenchmarkResult:
"""Confounder strengths calibrated against covariates that *were* observed.
Interpretation: ``cf_y`` and ``cf_d`` are the sensitivity parameters an unobserved
confounder would need in order to be "as important as" the benchmark covariates,
measured by how much dropping those covariates degrades the outcome regression and
the Riesz representer. ``delta_psi`` is how much the estimate actually moved when
they were dropped, and ``rho`` is the implied degree of adversity -- a value well
below 1 says the worst-case ``rho = 1`` is pessimistic for confounders like these.
Parameters
----------
estimand : str
Alias of the estimand benchmarked.
covariates : tuple of str
The observed covariates the strength is calibrated against.
cf_y : float
Gain in explained outcome variation from adjusting for those covariates,
``(R^2_long - R^2_short) / (1 - R^2_long)``, clipped to ``[0, 1]``. The raw value
is :attr:`gain_y`.
cf_d : float
Gain in the Riesz representer's second moment from adjusting for those
covariates, ``(nu2_long - nu2_short) / nu2_short``, clipped to ``[0, 1]``. The raw
value is :attr:`gain_d`. Both gains are the benchmark convention of Chernozhukov
et al. (arXiv:2112.13398v6, Appendix E.5, Remark 8) and of DoubleML. Remark 8
states that they serve as proxies for the shares that
:func:`omitted_variable_bounds` takes, so they enter it unchanged. Pass the
whole record as ``benchmark=``.
rho : float
The implied alignment of the two.
delta_psi : float
How far the estimate moved when they were dropped.
psi_long : float
Estimate with the benchmark covariates adjusted for.
psi_short : float
Estimate with them dropped.
sigma2_long : float
Residual outcome variance with them adjusted for.
sigma2_short : float
Residual outcome variance with them dropped.
nu2_long : float
Riesz second moment with them adjusted for.
nu2_short : float
Riesz second moment with them dropped.
random_state : int or None
Seed this benchmark ran under. Pass it back to :func:`benchmark` to obtain the
benchmark again. ``None`` only on a record built without a seed.
Attributes
----------
gain_y : float
gain_d : float
"""
estimand: str
covariates: tuple[str, ...]
cf_y: float
cf_d: float
rho: float
delta_psi: float
psi_long: float
psi_short: float
sigma2_long: float
sigma2_short: float
nu2_long: float
nu2_short: float
#: Seed the short refit ran under, resolved rather than requested: an explicit seed,
#: else the fit's own, else one drawn here. A benchmark refits, so an estimator
#: carrying no seed would otherwise give a different answer to the same question and
#: cache the first one.
random_state: int | None = None
@property
def gain_y(self) -> float:
"""The unclipped outcome gain, ``sigma2_short / sigma2_long - 1``.
A value above 1 says that a confounder as strong as the benchmark covariates would
explain all of the remaining outcome variation. ``cf_y`` then holds the largest
share, 1.
"""
if self.sigma2_long <= 0.0:
return float("inf")
return float(self.sigma2_short / self.sigma2_long - 1.0)
@property
def gain_d(self) -> float:
"""The unclipped Riesz gain, ``nu2_long / nu2_short - 1``.
A value of 1 or more gives ``cf_d = 1``, at which the bound is unbounded.
"""
if self.nu2_short <= 0.0:
return float("inf")
return float(self.nu2_long / self.nu2_short - 1.0)
[docs]
def to_dict(self) -> dict[str, Any]:
"""Return a JSON-compatible representation.
Returns
-------
dict
A JSON-compatible mapping of every reported field.
"""
return {
"estimand": self.estimand,
"covariates": ", ".join(self.covariates),
"cf_y": self.cf_y,
"cf_d": self.cf_d,
"gain_y": self.gain_y,
"gain_d": self.gain_d,
"rho": self.rho,
"delta_psi": self.delta_psi,
"psi_long": self.psi_long,
"psi_short": self.psi_short,
}
[docs]
def summary(self) -> str:
"""Return a printable summary.
Returns
-------
str
A printable report, one line per reported quantity.
"""
return "\n".join(
[
f"Benchmark for {self.estimand!r} against {list(self.covariates)}",
"-" * 48,
format_table(
["quantity", "with covariates", "without"],
[
["estimate", f"{self.psi_long:.5g}", f"{self.psi_short:.5g}"],
["sigma^2", f"{self.sigma2_long:.5g}", f"{self.sigma2_short:.5g}"],
["nu^2", f"{self.nu2_long:.5g}", f"{self.nu2_short:.5g}"],
],
),
"",
f"implied cf_y = {self.cf_y:.4f}, cf_d = {self.cf_d:.4f}, rho = {self.rho:.4f}",
*self._clip_lines(),
f"the estimate moved by {self.delta_psi:+.5g} when these covariates were dropped",
]
)
def _clip_lines(self) -> list[str]:
"""One line for each gain that the record clips at 1."""
lines = []
if self.gain_y > 1.0:
lines.append(
f"cf_y is the gain {self.gain_y:.4g} clipped at 1: a confounder this strong "
"explains all remaining outcome variation"
)
if self.gain_d >= 1.0:
lines.append(
f"cf_d is the gain {self.gain_d:.4g} clipped at 1: the bound is unbounded "
"for a confounder this strong"
)
return lines
def __repr__(self) -> str: # pragma: no cover - debugging aid
return self.summary()
def _benchmark_covariates(data: Any) -> tuple[str, ...]:
"""The covariates a benchmark can name, in fit order, with an encoded covariate once.
An encoded categorical covariate enters the fit as a block of indicator columns. One
indicator alone is not a covariate, so the column name stands for the whole block.
"""
sources = data.indicator_sources
return tuple(dict.fromkeys(sources.get(name, name) for name in data.covariate_names))
def _benchmark_columns(data: Any, names: tuple[str, ...]) -> tuple[str, ...]:
"""The design columns that the covariate ``names`` stand for, in request order."""
sources = data.indicator_sources
return tuple(
column
for name in names
for column in data.covariate_names
if sources.get(column, column) == name
)
def _benchmark_names(result: TMLEResult, covariates: Any) -> tuple[str, ...]:
"""The covariate names a ``covariates=`` value requests, one bare name included.
A malformed request is reported before any refusal: an empty one, an indicator column
of an encoded covariate, and a name the fit does not adjust for. :func:`benchmark`
and the ``benchmark`` method of :class:`~cleverly.assessment.SensitivityFacade` each
call this first, and :func:`benchmark_refusal` returns ``None`` where it raises.
Parameters
----------
result : TMLEResult
The fitted result whose covariates the names must be.
covariates : str or sequence of str
The requested covariate names. An encoded categorical covariate is named by its
column, not by one of its indicator columns.
Returns
-------
tuple of str
The requested names, in request order.
Raises
------
DataError
If the request names no covariate, names an indicator column, or names a
covariate the fit does not adjust for.
"""
names = tuple([covariates] if isinstance(covariates, str) else covariates)
if not names:
raise DataError(
"benchmark needs at least one covariate to drop; covariates= names none, and a "
"refit that drops nothing measures nothing"
)
sources = result.data.indicator_sources
for name in names:
if name in sources:
raise DataError(
f"{name!r} is an indicator column of the encoded covariate {sources[name]!r}. "
f"Name {sources[name]!r} to drop its whole encoded block"
)
fitted = _benchmark_covariates(result.data)
unknown = sorted(set(names).difference(fitted))
if unknown:
raise DataError(f"unknown covariates {unknown}; this fit adjusts for {list(fitted)}")
return names
def benchmark_refusal(result: TMLEResult, covariates: Any = None) -> Refusal | None:
"""Return the refusal a :func:`benchmark` request meets for its ``covariates``.
:func:`benchmark` raises the sentence this function returns, and the ``benchmark``
capability row of :class:`~cleverly.assessment.SensitivityFacade` quotes it for the
same request. The short model is a refit without the named covariates. A request
that names every covariate of the fit leaves the refit nothing to adjust for, and
this package fits no TMLE without a covariate, so the request is refused. The
unadjusted comparison is well posed, and no estimator here computes it.
Parameters
----------
result : TMLEResult
A fitted point-treatment result.
covariates : str, sequence of str, or None
The requested covariate names. ``None`` stands for the omitted argument, which
no value can lift on a fit with one covariate.
Returns
-------
Refusal or None
The refusal, or ``None`` when the covariates leave one to adjust for or
form a malformed request. :func:`benchmark` reports a malformed request itself,
before any refusal.
"""
if result.data.known_treatment is not None:
return Refusal(
RefusalReason.INVALID_DECLARATION,
"benchmark compares the treatment mechanism with and without the named "
"covariates, and this fit's data declares the mechanism known, so dropping a "
"covariate cannot change it. A known design mechanism leaves no omitted "
"confounder of the treatment to calibrate against.",
)
fitted = _benchmark_covariates(result.data)
try:
dropped = set(fitted[:1] if covariates is None else _benchmark_names(result, covariates))
except DataError:
return None
if any(name not in dropped for name in fitted):
return None
return Refusal(
RefusalReason.IMPLEMENTATION_GAP,
f"benchmark cannot drop every covariate of this fit, {list(fitted)}. The short "
"model would adjust for nothing, and a fit needs at least one covariate.",
alternative=(
"Name a proper subset of the covariates, which exists only when the fit "
"adjusts for two or more."
),
roadmap="X65",
)
[docs]
def benchmark(
result: TMLEResult,
covariates: Any,
*,
estimand: str = "ate",
nu2_estimator: str = "auto",
random_state: int | None = None,
) -> BenchmarkResult:
"""Calibrate ``cf_y`` and ``cf_d`` against observed covariates.
Refits the whole model *without* the named covariates (the "short" model) and
compares it with the full fit. The resulting gain statistics say how strong a
confounder like the dropped ones would be, on the ``cf_y``/``cf_d`` scale, which
turns an abstract sensitivity parameter into a concrete comparison.
Note this is a genuine refit, so it costs about as much as the original fit.
Parameters
----------
result : TMLEResult
A fitted result.
covariates : sequence of str
Observed covariates to calibrate against. An encoded categorical covariate is
named by its column, and the refit drops its whole block of indicator columns.
estimand : str
Alias to benchmark.
nu2_estimator : {"auto", "doubly_robust", "plugin"}
Which estimator of the Riesz second moment to use.
random_state : int or None
Seed for the short refit. ``None`` uses the seed the fit was run with. A fit
with no seed draws one. Either way the result records it under ``random_state``,
and passing that value back repeats the benchmark.
Returns
-------
BenchmarkResult
The sensitivity parameters a confounder as important as those covariates
would need, and how far the estimate moved without them.
Raises
------
DataError
If ``covariates`` is empty, names an indicator column of an encoded covariate, or
names a covariate the fit does not adjust for. The argument is checked before any
refusal.
ValueError
If ``nu2_estimator`` is not one of :data:`NU2_ESTIMATORS`.
CapabilityError
If the result carries no fitted estimator, on every refusal
:func:`sensitivity_elements` raises for the full fit or the short refit, and on
the refusal of ``benchmark_refusal``. The full fit and the covariates are
checked before the refit runs.
"""
names = _benchmark_names(result, covariates)
estimator = result.estimator
if estimator is None:
raise CapabilityError(
"benchmark needs the fitted estimator that produced the result",
reason=RefusalReason.INVALID_DECLARATION,
)
# The short model is a refit, so it carries the same reproducibility question a
# refutation does: an estimator with no ``random_state`` redraws its folds every time,
# and the result is cached on the fit and survives ``save``. Resolve one seed, run the
# refit under it, and report it. Same convention as ``cleverly.validation.refute``.
seed = resolve_assessment_seed(result, random_state)
long_elements = sensitivity_elements(result, estimand, nu2_estimator=nu2_estimator)
refusal = benchmark_refusal(result, names)
if refusal is not None:
raise CapabilityError.from_refusal(refusal)
short_data = result.data.without_covariates(_benchmark_columns(result.data, names))
short_result = estimator.refit(
short_data, intermediate_value=result.intermediate_value, random_state=seed
)
short_elements = sensitivity_elements(short_result, estimand, nu2_estimator=nu2_estimator)
var_y = float(np.var(result.data.outcome[result.data.observed]))
r2_long = 1.0 - long_elements.sigma2 / var_y
r2_short = 1.0 - short_elements.sigma2 / var_y
r2_riesz = short_elements.nu2 / long_elements.nu2
cf_y = float(np.clip((r2_long - r2_short) / (1.0 - r2_long), 0.0, 1.0))
cf_d = float(np.clip((1.0 - r2_riesz) / r2_riesz, 0.0, 1.0)) if r2_riesz > 0 else 1.0
delta = short_result[estimand].psi - result[estimand].psi
var_g = short_elements.sigma2 - long_elements.sigma2
var_riesz = long_elements.nu2 - short_elements.nu2
if var_g > 0 and var_riesz > 0:
rho = float(np.clip(abs(delta) / np.sqrt(var_g * var_riesz), 0.0, 1.0)) * float(
np.sign(delta)
)
else:
rho = float(np.sign(delta))
return BenchmarkResult(
estimand=estimand,
covariates=names,
cf_y=cf_y,
cf_d=cf_d,
rho=rho,
delta_psi=delta,
psi_long=result[estimand].psi,
psi_short=short_result[estimand].psi,
sigma2_long=long_elements.sigma2,
sigma2_short=short_elements.sigma2,
nu2_long=long_elements.nu2,
nu2_short=short_elements.nu2,
random_state=seed,
)
[docs]
def robustness_value(
result: TMLEResult,
estimand: str = "ate",
*,
rho: float = 1.0,
level: float = 0.95,
null_hypothesis: float = 0.0,
nu2_estimator: str = "auto",
) -> dict[str, Any]:
"""The confounding strength that would explain the effect away.
Returns ``{"rv": ..., "rva": ...}``: the value of ``cf_y = cf_d`` at which the point
estimate reaches ``null_hypothesis``, and the value at which its confidence bound
does. ``None`` means no equal strength below one reaches the null. A confidence
limit that already crosses the null has ``rva = 0``. This is the single most useful
number in this module, because it requires no
guess about how strong an unmeasured confounder might be -- it reports the
threshold and lets the reader judge whether it is plausible.
The confidence bound is read off the estimate's influence curve, so at a
non-inferential status ``"rva"`` takes the name
:func:`~cleverly.inference.influence.spread_name` gives it, ``"rv_plugin_interval"``,
and the mapping opens with ``"inference"`` naming the status, as
:meth:`SensitivityBounds.to_dict` does. Under ``nu2_estimator="plugin"`` no
derivation gives the confidence bound's standard error, so the mapping carries
``"nu2_estimator"`` and no confidence-limit value under any name (docs/roadmap.md F26).
Parameters
----------
result : TMLEResult
A fitted result.
estimand : str
Alias to report on.
rho : float
How adversarially the two sensitivity parameters are aligned.
level : float
Coverage level used for the confidence-limit value.
null_hypothesis : float
Value the strength is measured against.
nu2_estimator : {"auto", "doubly_robust", "plugin"}
Which estimator of the Riesz second moment to use.
Returns
-------
dict of str to Any
The strength that moves the point estimate to the null, the one that moves the
confidence limit there, and ``max_bias``. A threshold is ``None`` when unreachable.
At a non-inferential status, also the
status under ``"inference"``. Under the plug-in estimator, ``"nu2_estimator"``
in place of the confidence-limit value.
Raises
------
ValueError
If ``nu2_estimator`` is not one of :data:`NU2_ESTIMATORS`, or if ``|rho|``
exceeds one, or if ``level`` lies outside ``[0.5, 1)``.
CapabilityError
On every refusal :func:`sensitivity_elements` raises: a fit the rule table
refuses, an estimand the bound does not cover, or a doubly robust estimate of
``nu^2`` that is not positive.
"""
elements = sensitivity_elements(result, estimand, nu2_estimator=nu2_estimator)
estimate = result[estimand]
rv, rva = _robustness_values(
result, elements, estimate.psi, rho, level, null_hypothesis, estimate.reference_df
)
status = estimate.inference
values: dict[str, Any] = {"inference": status} if is_labelled(status) else {}
if elements.psi_max_bias is None:
values |= {
"nu2_estimator": elements.nu2_estimator,
"rv": rv,
"max_bias": elements.max_bias,
}
if elements.max_bias == 0.0 and elements.nu2_estimator == "doubly_robust":
values["limit_refusal"] = _ZERO_BIAS_LIMITS_REFUSAL
return values
values |= {"rv": rv, spread_name("rva", status): rva, "max_bias": elements.max_bias}
return values
def contour_data(
result: TMLEResult,
estimand: str = "ate",
*,
rho: float = 1.0,
grid_size: int = 20,
grid_bounds: tuple[float, float] = (0.15, 0.15),
bound: str = "lower",
nu2_estimator: str = "auto",
) -> Any:
"""A ``cf_d`` x ``cf_y`` grid of bias-adjusted bounds, for a contour plot.
Returned as a long-format frame (``cf_d``, ``cf_y``, ``value``) rather than a plot, so
it can be rendered with whatever plotting stack the caller already uses.
Parameters
----------
result : TMLEResult
A fitted result.
estimand : str
Alias to build the grid for.
rho : float
How adversarially the two sensitivity parameters are aligned.
grid_size : int
Points along each axis of the grid.
grid_bounds : tuple of float
Largest ``cf_d`` and largest ``cf_y`` the grid reaches.
bound : {"lower", "upper"}
Which end of the bias-adjusted interval each cell reports.
nu2_estimator : {"auto", "doubly_robust", "plugin"}
Which estimator of the Riesz second moment to use.
Returns
-------
DataFrame
One row per grid cell, in the frame library the fit was given, with columns
``cf_d``, ``cf_y``, and ``value``.
Raises
------
ValueError
If ``bound`` is not ``"lower"`` or ``"upper"``, if ``nu2_estimator`` is not one
of :data:`NU2_ESTIMATORS`, if a ``cf_y`` value lies outside ``[0, 1]`` or a
``cf_d`` value outside ``[0, 1)``, or if
``|rho|`` exceeds one.
CapabilityError
On every refusal :func:`sensitivity_elements` raises: a fit the rule table
refuses, an estimand the bound does not cover, or a doubly robust estimate of
``nu^2`` that is not positive.
"""
if bound not in ("lower", "upper"):
raise ValueError(f"bound must be 'lower' or 'upper'; got {bound!r}")
elements = sensitivity_elements(result, estimand, nu2_estimator=nu2_estimator)
psi = result[estimand].psi
sign = -1.0 if bound == "lower" else 1.0
cf_d_grid = np.linspace(0.0, grid_bounds[0], grid_size)
cf_y_grid = np.linspace(0.0, grid_bounds[1], grid_size)
rows_d, rows_y, values = [], [], []
for cf_d in cf_d_grid:
for cf_y in cf_y_grid:
strength = _confounding_strength(float(cf_y), float(cf_d), rho)
rows_d.append(float(cf_d))
rows_y.append(float(cf_y))
values.append(psi + sign * strength * elements.max_bias)
return result.data.frame_like({"cf_d": rows_d, "cf_y": rows_y, "value": values})