r"""Refutation tests: try to break the fit on purpose.
Diagnostics tell you about the fit you have. Refutation tests -- in the spirit of
DoWhy's refuters -- perturb the *problem* in ways whose correct answer is known in
advance, and check that the estimator gives that answer. They cannot prove the
identification assumptions hold, but they do catch real defects: leakage, a
mis-specified estimand, an implementation that is quietly picking up an artefact.
The tests, and what a failure means:
``placebo``
Replace treatment with a random permutation of itself, destroying any real
effect while preserving its marginal distribution. The estimate should be
indistinguishable from the parameter's no-effect value. A non-null placebo
estimate can mean leakage between folds or a misdefined estimand.
``random_common_cause``
Add an irrelevant random covariate. The estimate should barely move; the
adjustment set is supposed to be robust to noise. A large shift means the
nuisance models are unstable at this sample size.
``subset``
Refit on random subsamples. The estimates should scatter around the full-sample
estimate by roughly the standard error. Scatter far exceeding it points to
influential observations -- usually the same handful of units the positivity
diagnostics flag.
``negative_control_outcome``
Refit with an outcome that the treatment cannot plausibly affect. Under a valid,
comparable negative-control design, a non-null estimate can reveal residual bias.
It can also reflect a bad control or a bias that does not affect the primary outcome.
A null estimate does not prove that unmeasured confounding is absent.
``dummy_outcome``
Generate Gaussian noise independent of treatment and adjustment variables, then
fully refit. Its declared effect is zero; exclusion from the empirical refit
distribution reveals a pipeline that manufactures signal from a null outcome.
``simulated_outcome``
Generate a Gaussian outcome from a recorded adjustment function, additive treatment
term, and noise law, then fully refit. The declared treatment effect must lie inside
the empirical refit distribution. Failure reveals that the pipeline cannot recover
a known effect under its first registered generated-outcome process.
``bootstrap_measurement_error``
Draw a bootstrap sample, add declared measurement error to the selected adjustment
variables, then fully refit. Mild measurement error should move the estimate by less
than sampling noise, so the original estimate must lie inside the empirical refit
distribution. Exclusion means the fit rests on adjustment variables it cannot
tolerate any error in. A draw that fails to sample, perturb, or refit is retained as
a failure rather than raised, because a bootstrap draw fails per sample.
Each test refits the model, so a full run costs several times a single fit.
"""
from __future__ import annotations
from collections.abc import Callable, Sequence
from dataclasses import dataclass, field, replace
from typing import TYPE_CHECKING, Any, Literal
import numpy as np
from .._inference_status import InferenceStatus
from ..exceptions import NO_ALTERNATIVE, CapabilityError, Refusal, RefusalReason
from ..inference.bootstrap import Resampling, _bootstrap_design
from ..inference.influence import spread_name
from ..utils.frames import emit_frame
from ..utils.random import resolve_assessment_seed
from ..utils.text import format_table
from .simulation import ReplicationFailure
if TYPE_CHECKING: # pragma: no cover - typing only
from ..estimators.base import TMLEResult
__all__ = [
"DEFAULT_OUTCOME_REPLICATES",
"DEFAULT_TESTS",
"BootstrapMeasurementError",
"EmpiricalInclusionRule",
"EmpiricalRefitRecord",
"GaussianAdjustmentOutcome",
"GaussianIndependentOutcome",
"GaussianNoise",
"RefutationResult",
"RefutationTest",
"RelativeGaussianNoise",
"refute",
]
#: Tests run by default. ``negative_control_outcome`` is excluded because it needs a
#: control outcome only the analyst can supply.
DEFAULT_TESTS: tuple[str, ...] = ("placebo", "random_common_cause", "subset")
#: Generated outcomes and measurement-error refits need an empirical distribution.
#: The other perturbation tests use five refits by default.
DEFAULT_OUTCOME_REPLICATES = 100
_DEFAULT_PERTURBATION_REPLICATES = 5
[docs]
@dataclass(frozen=True)
class GaussianNoise:
"""Declare one Gaussian noise law for a generated outcome.
Parameters
----------
mean : float
Mean of each independent noise draw.
standard_deviation : float
Positive standard deviation of each independent noise draw.
"""
mean: float = 0.0
standard_deviation: float = 1.0
def __post_init__(self) -> None:
if not np.isfinite(self.mean):
raise ValueError("Gaussian noise mean must be finite")
if not np.isfinite(self.standard_deviation) or self.standard_deviation <= 0.0:
raise ValueError("Gaussian noise standard_deviation must be finite and positive")
[docs]
def draw(self, rng: np.random.Generator, size: int) -> np.ndarray:
"""Draw independent values from the declared law.
Parameters
----------
rng : numpy.random.Generator
Generator that supplies the draw.
size : int
Number of values to draw.
Returns
-------
ndarray
Independent Gaussian values.
"""
return np.asarray(rng.normal(self.mean, self.standard_deviation, size=size), dtype=float)
[docs]
@dataclass(frozen=True)
class RelativeGaussianNoise:
"""Declare mean-zero numeric noise relative to a sampled variable.
Parameters
----------
standard_deviation : float
Nonnegative multiple of the variable's bootstrap-sample standard deviation.
"""
standard_deviation: float = 0.1
def __post_init__(self) -> None:
if not np.isfinite(self.standard_deviation) or self.standard_deviation < 0.0:
raise ValueError("relative Gaussian standard_deviation must be finite and nonnegative")
[docs]
def draw(self, rng: np.random.Generator, values: Any) -> np.ndarray:
"""Draw noise on the declared bootstrap-sample scale.
Parameters
----------
rng : numpy.random.Generator
Generator that supplies the draw.
values : array-like
Sampled numeric variable whose population standard deviation sets the scale.
Returns
-------
ndarray
One mean-zero Gaussian draw per value.
"""
sample = np.asarray(values, dtype=float).reshape(-1)
scale = self.standard_deviation * float(np.std(sample, ddof=0))
return np.asarray(rng.normal(0.0, scale, size=sample.size), dtype=float)
[docs]
@dataclass(frozen=True)
class BootstrapMeasurementError:
"""Declare bootstrap measurement-error validation for adjustment variables.
Parameters
----------
variables : tuple of str
Original logical adjustment-variable names to perturb.
numeric_noise : RelativeGaussianNoise
Relative law for each selected numeric variable.
categorical_change_probability : float
Probability that a selected categorical value changes to another level.
resampling : {"auto", "iid", "cluster"}
Bootstrap unit. ``"auto"`` uses clusters when the fit declares them.
"""
variables: tuple[str, ...]
numeric_noise: RelativeGaussianNoise = RelativeGaussianNoise()
categorical_change_probability: float = 0.1
resampling: Resampling = "auto"
def __post_init__(self) -> None:
if isinstance(self.variables, str):
raise ValueError("bootstrap measurement error variables must be a sequence of names")
try:
selected = tuple(self.variables)
except TypeError:
# ``tuple`` on a non-iterable raises ``TypeError``, which would escape the
# module's ``ValueError`` convention for a rejected declaration.
raise ValueError(
"bootstrap measurement error variables must be a sequence of names"
) from None
object.__setattr__(self, "variables", selected)
if not self.variables:
raise ValueError("bootstrap measurement error variables must not be empty")
if len(set(self.variables)) != len(self.variables):
raise ValueError("bootstrap measurement error variables must not repeat")
if not all(isinstance(name, str) and name for name in self.variables):
raise ValueError("bootstrap measurement error variables must be nonempty names")
if type(self.numeric_noise) is not RelativeGaussianNoise:
raise ValueError("numeric_noise must be a RelativeGaussianNoise declaration")
probability = self.categorical_change_probability
if not np.isfinite(probability) or not 0.0 <= probability <= 1.0:
raise ValueError("categorical_change_probability must be between zero and one")
if self.resampling not in ("auto", "iid", "cluster"):
raise ValueError("resampling must be 'auto', 'iid', or 'cluster'")
[docs]
@dataclass(frozen=True)
class GaussianIndependentOutcome:
"""Declare an outcome independent of treatment and adjustment variables.
Parameters
----------
noise : GaussianNoise
Gaussian law that generates every outcome value.
"""
noise: GaussianNoise = GaussianNoise()
name: str = field(default="gaussian_independent_noise", init=False)
adjustment_function: str = field(default="none", init=False)
treatment_term: str = field(default="none", init=False)
family: str = field(default="gaussian", init=False)
known_effect: float = field(default=0.0, init=False)
@property
def noise_law(self) -> GaussianNoise:
"""Return the recorded Gaussian noise declaration."""
return self.noise
[docs]
def draw(self, data: Any, *, seed: int) -> np.ndarray:
"""Draw one independent outcome vector.
Parameters
----------
data : CausalData
Analysis rows that determine the output length.
seed : int
Child seed for this draw.
Returns
-------
ndarray
One generated outcome per analysis row.
"""
return self.noise.draw(np.random.default_rng(seed), data.n)
[docs]
@dataclass(frozen=True)
class GaussianAdjustmentOutcome:
r"""Declare ``f(W) + effect * A + epsilon`` with Gaussian noise.
The adjustment function standardizes each nonconstant covariate, sums the standardized
columns, divides by the square root of how many columns it summed, and multiplies the
result by ``adjustment_scale``. The root-count divisor holds the adjustment standard
deviation at about ``adjustment_scale`` for any number of approximately independent
covariates, so the confounding signal does not shrink against the fixed noise scale as
the adjustment set widens. Correlated columns raise that standard deviation instead,
which only makes a failure to adjust easier to detect. The known binary arm contrast is
``effect`` because the treatment term is additive.
Parameters
----------
effect : float
Declared effect of treatment code one against code zero.
adjustment_scale : float
Multiplier on the root-count normalized standardized covariate sum.
noise : GaussianNoise
Gaussian law for ``epsilon``.
"""
effect: float = 1.0
adjustment_scale: float = 1.0
noise: GaussianNoise = GaussianNoise()
name: str = field(default="gaussian_adjustment_dependent", init=False)
adjustment_function: str = field(
default="scaled_standardized_covariate_sum_over_root_count", init=False
)
treatment_term: str = field(default="effect_times_treatment_code", init=False)
family: str = field(default="gaussian", init=False)
def __post_init__(self) -> None:
if not np.isfinite(self.effect):
raise ValueError("Gaussian adjustment outcome effect must be finite")
if not np.isfinite(self.adjustment_scale):
raise ValueError("Gaussian adjustment outcome adjustment_scale must be finite")
@property
def known_effect(self) -> float:
"""Return the derived code-one against code-zero mean contrast."""
return float(self.effect)
@property
def noise_law(self) -> GaussianNoise:
"""Return the recorded Gaussian noise declaration."""
return self.noise
[docs]
def adjustment(self, covariates: Any) -> np.ndarray:
"""Evaluate the declared adjustment function.
Parameters
----------
covariates : array-like
Covariate matrix with one row per observation.
Returns
-------
ndarray
Adjustment value for each row.
"""
values = np.asarray(covariates, dtype=float)
center = np.mean(values, axis=0)
scale = np.std(values, axis=0)
active = scale > 0.0
if not np.any(active):
return np.zeros(values.shape[0], dtype=float)
standardized = (values[:, active] - center[active]) / scale[active]
# Sum over root count rather than mean: a mean of p standardized columns has
# standard deviation about 1/sqrt(p), so the confounding signal would vanish
# against the fixed noise scale on a wide adjustment set.
divisor = np.sqrt(float(np.count_nonzero(active)))
return np.asarray(
self.adjustment_scale * np.sum(standardized, axis=1) / divisor, dtype=float
)
[docs]
def draw(self, data: Any, *, seed: int) -> np.ndarray:
"""Draw one adjustment-dependent outcome vector.
Parameters
----------
data : CausalData
Analysis rows that supply covariates and treatment.
seed : int
Child seed for the Gaussian noise.
Returns
-------
ndarray
One generated outcome per analysis row.
"""
noise = self.noise.draw(np.random.default_rng(seed), data.n)
# A row with an unrecorded treatment has no code to add an effect to. Its
# generated outcome is never read, because the composite indicator is zero there,
# so the arm-0 value keeps the vector finite and leaves the known effect intact.
treatment = np.where(np.isfinite(data.treatment), data.treatment, 0.0)
return np.asarray(
self.adjustment(data.covariates) + self.effect * treatment + noise,
dtype=float,
)
[docs]
@dataclass(frozen=True)
class EmpiricalInclusionRule:
"""Declare a two-sided empirical inclusion rule for generated outcomes.
Parameters
----------
alpha : float
Two-sided tail probability.
minimum_draws : int
Minimum successful draw count needed to apply the empirical rule. The rule refuses
a count below ``2 / alpha``, where the declared level is unreachable.
failure_policy : {"fail"}
Policy for a failed refit. The first catalog always fails if any refit fails.
"""
alpha: float = 0.05
minimum_draws: int = 40
failure_policy: Literal["fail"] = "fail"
def __post_init__(self) -> None:
if not 0.0 < self.alpha < 1.0:
raise ValueError("empirical inclusion alpha must be between zero and one")
if self.minimum_draws < 1:
raise ValueError("empirical inclusion minimum_draws must be positive")
# The smallest positive value ``pvalue`` can return on n draws is 2/n. When that
# exceeds alpha the rule can only reject a sample whose draws all fall on one side
# of the declared effect, which leaves it with almost no power: it stops reading the
# tail of the refit distribution and reads only whether the truth is outside its
# whole range.
needed = int(np.ceil(2.0 / self.alpha - 1e-9))
if self.minimum_draws < needed:
raise ValueError(
f"empirical inclusion needs minimum_draws >= {needed} at alpha={self.alpha:g}; "
f"{self.minimum_draws} draws cannot reach a two-sided p-value at or below "
"alpha except on an all-on-one-side sample"
)
if self.failure_policy != "fail":
raise ValueError("generated-outcome refutations support only failure_policy='fail'")
[docs]
def pvalue(self, values: Sequence[float], truth: float) -> float:
"""Return a two-sided empirical rank p-value with inclusive half-ties.
Parameters
----------
values : sequence of float
Successful refitted estimates.
truth : float
Declared process effect.
Returns
-------
float
Two-sided empirical rank probability. An empty sample returns ``0.0``, which is
a rank over nothing rather than a rejection. :meth:`evaluate` reports ``None``
for that case, and it is the value a report carries.
"""
estimates = np.asarray(values, dtype=float)
if estimates.size == 0:
return 0.0
below = float(np.count_nonzero(estimates < truth))
above = float(np.count_nonzero(estimates > truth))
ties = float(estimates.size - below - above)
return float(min(1.0, 2.0 * min(below + 0.5 * ties, above + 0.5 * ties) / estimates.size))
@property
def minimum_draw_count(self) -> int:
"""Return the minimum successful draw count under this rule."""
return self.minimum_draws
[docs]
def evaluate(
self,
values: Sequence[float],
truth: float,
failures: Sequence[ReplicationFailure],
) -> tuple[bool, float | None, str]:
"""Apply the declared draw-count, failure, and empirical-tail gates.
Parameters
----------
values : sequence of float
Successful refitted estimates.
truth : float
Declared process effect.
failures : sequence of ReplicationFailure
Failed refits retained by the operation.
Returns
-------
tuple of bool, float or None, str
Pass status, empirical probability, and readable detail. The probability is
``None`` when no draw succeeded, because an empty sample carries no rank
information and a reported zero would read as a strong rejection.
"""
# ``None`` rather than a number the empty sample cannot support.
probability = self.pvalue(values, truth) if len(values) else None
if failures:
return (
False,
probability,
f"{len(failures)} refit(s) failed; failure_policy='fail' keeps every failure",
)
if len(values) < self.minimum_draws:
return (
False,
probability,
f"only {len(values)} successful draw(s); the rule requires {self.minimum_draws}",
)
probability = self.pvalue(values, truth)
return (
probability > self.alpha,
probability,
f"two-sided empirical p={probability:.5g} with inclusive half-ties "
f"(alpha={self.alpha:g})",
)
[docs]
@dataclass(frozen=True)
class EmpiricalRefitRecord:
"""One successful empirical refit.
Parameters
----------
replicate : int
Zero-based draw index.
seed : int
Child seed for this draw. It seeds the estimator refit, and the perturbation the
operation applies before it: outcome generation, or bootstrap measurement error.
estimand : str
Structured result alias selected for the refit.
estimate : float
Refitted point estimate.
std_error : float
Refitted standard error, read through
:attr:`~cleverly.ParameterEstimate.plugin_std_error`. On a fit that supplies no
inference it is that diagnostic, and not a standard error.
family : str
Outcome family supplied by the refitted result.
"""
replicate: int
seed: int
estimand: str
estimate: float
std_error: float
family: str
def _format_number(value: float | None) -> str:
"""Format one optional number for a summary line.
Parameters
----------
value : float or None
Number to format.
Returns
-------
str
Five significant figures, and ``"-"`` for a missing value.
"""
return "-" if value is None else f"{value:.5g}"
[docs]
@dataclass(frozen=True)
class RefutationTest:
"""The outcome of one refutation test.
Parameters
----------
name : str
Which refutation was run.
estimand : str
Alias the test was run for.
original : float
The estimate before the refutation.
values : tuple of float
Estimates the refutation produced.
expectation : str
What those values should look like if the fit is sound.
passed : bool
Whether they did.
detail : str
What was seen, in the test's own terms.
standard_errors : tuple of float
Standard errors from successful empirical refits.
declaration : object or None
Immutable declaration, for an empirical refit test. It holds the outcome process,
or the measurement-error request.
rule : EmpiricalInclusionRule or None
Recorded inclusion rule, for an empirical refit test.
records : tuple of EmpiricalRefitRecord
Successful empirical refit draws.
failures : tuple of ReplicationFailure
Failed empirical draws retained by the shared failure contract. Every operation
retains a refit failure. ``bootstrap_measurement_error`` also retains a draw it
cannot sample or perturb.
family : str or None
Family supplied by successful refitted results.
empirical_pvalue : float or None
Two-sided empirical rank p-value for the declared effect. ``None`` when no draw
succeeded, and on an older refuter.
declared_effect : float or None
Effect the declared process implies for this test's own contrast direction, and
``None`` for an older refuter. It carries the sign of the parameter key: a fit
that reports code zero against code one declares the negated process effect.
requested_draws : int or None
Requested empirical draw count.
resampling : str or None
Resolved bootstrap mode for a bootstrap-based refuter.
inference : str
Inference status of the refitted estimates. One of
:data:`~cleverly.inference.influence.InferenceStatus`, which the
:doc:`inference reference </technical-reference/inference>` lists. On a fit
that supplies no inference, ``standard_errors`` and each record's ``std_error``
hold the plug-in diagnostic
:attr:`~cleverly.ParameterEstimate.plugin_std_error`, and :meth:`to_frame`
publishes that column as ``plugin_std_error``.
mean_kind : str
``"geometric"`` for a ratio placebo or noise refutation tested on the log
scale. Other tests use the arithmetic mean. The underlying ``values`` remain
on the reported scale.
spread_kind : str
``"log"`` for ratio subset scatter tested on the log scale. Other reports use
reported-scale scatter. The underlying ``values`` remain on the reported scale.
"""
name: str
estimand: str
original: float
values: tuple[float, ...]
expectation: str
passed: bool
detail: str
standard_errors: tuple[float, ...] = ()
declaration: Any = None
rule: EmpiricalInclusionRule | None = None
records: tuple[EmpiricalRefitRecord, ...] = ()
failures: tuple[ReplicationFailure, ...] = ()
family: str | None = None
empirical_pvalue: float | None = None
declared_effect: float | None = None
requested_draws: int | None = None
resampling: str | None = None
inference: InferenceStatus = "influence_curve"
mean_kind: Literal["arithmetic", "geometric"] = "arithmetic"
spread_kind: Literal["reported", "log"] = "reported"
@property
def mean(self) -> float:
"""Return the displayed mean of refitted estimates.
Returns
-------
float
Arithmetic mean by default, or geometric mean for a ratio placebo or noise
refutation.
``nan`` when no replicate is finite. Failed replicates are dropped.
"""
finite = np.asarray([v for v in self.values if np.isfinite(v)])
if not finite.size:
return float("nan")
if self.mean_kind == "geometric":
return float(np.exp(np.mean(np.log(finite))))
return float(finite.mean())
@property
def spread(self) -> float:
"""Return how far this test's replicates spread on the verdict scale.
Returns
-------
float
Sample standard deviation over the finite replicates, on the log scale
for a ratio subset test. Returns ``nan`` with fewer than two values.
"""
finite = np.asarray([v for v in self.values if np.isfinite(v)])
if self.spread_kind == "log":
finite = np.log(finite)
return float(finite.std(ddof=1)) if finite.size > 1 else float("nan")
@property
def n_failed(self) -> int:
"""Return the number of retained failed empirical draws."""
return len(self.failures)
@property
def n_replicates(self) -> int:
"""Return the number of requested draws represented by this test."""
if self.declaration is not None:
return len(self.records) + len(self.failures)
return len(self.values)
@property
def process(self) -> Any:
"""Return the empirical refit declaration, or ``None`` for older refuters."""
return self.declaration
@property
def child_seeds(self) -> tuple[int, ...]:
"""Return every successful and failed draw seed in draw-index order."""
indexed = [(item.replicate, item.seed) for item in self.records]
indexed.extend((item.replicate, item.seed) for item in self.failures)
return tuple(seed for _, seed in sorted(indexed))
[docs]
def to_frame(self, data: Any = None, *, backend: str | None = None) -> Any:
"""Return one row per empirical refit draw.
Parameters
----------
data : Any
Dataframe or fitted container whose backend to match.
backend : str or None
Backend name used when ``data`` is omitted.
Returns
-------
dataframe
Successful estimates and standard errors alongside retained failures. The
standard-error column is ``plugin_std_error`` when :attr:`inference` is a
diagnostic status.
"""
successes = {item.replicate: item for item in self.records}
failures = {item.replicate: item for item in self.failures}
indices = sorted((*successes, *failures))
error_column = spread_name("std_error", self.inference)
payload: dict[str, list[Any]] = {
"test": [],
"replicate": [],
"seed": [],
"estimate": [],
error_column: [],
"family": [],
"error_type": [],
"message": [],
}
for index in indices:
record = successes.get(index)
failure = failures.get(index)
# Every index comes from one of the two maps, so exactly one lookup can miss.
# Read the seed off whichever hit rather than narrowing with an ``assert``,
# which ``-O`` strips.
seed = record.seed if record is not None else failures[index].seed
payload["test"].append(self.name)
payload["replicate"].append(index)
payload["seed"].append(seed)
payload["estimate"].append(None if record is None else record.estimate)
payload[error_column].append(None if record is None else record.std_error)
payload["family"].append(None if record is None else record.family)
payload["error_type"].append(None if failure is None else failure.error_type)
payload["message"].append(None if failure is None else failure.message)
return emit_frame(payload, data, backend=backend)
[docs]
@dataclass(frozen=True)
class RefutationResult:
"""All refutation tests run on a fit.
Parameters
----------
tests : tuple of RefutationTest
One record per refutation run.
estimand : str
Alias the tests were run for.
backend : str or None
Dataframe backend :meth:`to_frame` returns when ``data`` is omitted.
random_state : int or None
Seed this report ran under. Pass it back to :func:`refute` to obtain the report
again, whether or not the fit carries a seed of its own. ``None`` only on a
report built by hand.
"""
tests: tuple[RefutationTest, ...]
estimand: str
#: Name of the dataframe backend the fit's data arrived in, so that
#: :meth:`to_frame` honours "results come back in the backend you passed in"
#: without a caller having to thread the container back in by hand.
backend: str | None = None
#: Seed this report ran under, resolved rather than requested: an explicit seed, else
#: the fit's own, else one drawn here. It pins the perturbations and the refits they
#: feed, so the report is enough to obtain the report again.
random_state: int | None = None
@property
def passed(self) -> bool:
"""Return whether every required check passed."""
return all(test.passed for test in self.tests)
def __bool__(self) -> bool:
return self.passed
def __getitem__(self, name: str) -> RefutationTest:
for test in self.tests:
if test.name == name:
return test
raise KeyError(f"no test named {name!r}; have {[t.name for t in self.tests]}")
[docs]
def to_frame(self, data: Any = None) -> Any:
"""Return tabular output in the input dataframe backend.
Parameters
----------
data : Any
A dataframe or fitted container whose backend to match. ``None`` uses the
backend recorded on this object.
Returns
-------
dataframe
One row per refutation test.
"""
payload = {
"test": [test.name for test in self.tests],
"estimand": [test.estimand for test in self.tests],
"original": [test.original for test in self.tests],
"refuted_mean": [test.mean for test in self.tests],
"refuted_sd": [test.spread for test in self.tests],
"declared_effect": [test.declared_effect for test in self.tests],
"family": [test.family for test in self.tests],
"successful_draws": [len(test.records) for test in self.tests],
"failed_draws": [test.n_failed for test in self.tests],
"empirical_pvalue": [test.empirical_pvalue for test in self.tests],
"requested_draws": [test.requested_draws for test in self.tests],
"resampling": [test.resampling for test in self.tests],
"expectation": [test.expectation for test in self.tests],
"passed": [test.passed for test in self.tests],
}
return emit_frame(payload, data, backend=self.backend)
[docs]
def draws_frame(self, name: str, data: Any = None) -> Any:
"""Return successful and failed draws for one empirical refit test.
Parameters
----------
name : str
Name of a generated-outcome or measurement-error refutation in this report.
data : Any
Dataframe or fitted container whose backend to match. ``None`` uses the
backend recorded by the report.
Returns
-------
dataframe
One row per requested empirical draw.
"""
test = self[name]
if test.declaration is None:
raise ValueError(f"refutation {name!r} has no empirical draw records")
return test.to_frame(data, backend=self.backend)
[docs]
def summary(self) -> str:
"""Return a printable summary.
Returns
-------
str
A printable table, one line per refutation test.
"""
lines = [
f"Refutation tests for {self.estimand!r}",
"-" * 34,
format_table(
["test", "original", "refuted (mean)", "sd", "expected", "ok"],
[
[
test.name,
f"{test.original:.5g}",
f"{test.mean:.5g}",
f"{test.spread:.4g}" if np.isfinite(test.spread) else "-",
test.expectation,
"yes" if test.passed else "NO",
]
for test in self.tests
],
),
]
for test in self.tests:
if test.declaration is not None:
lines.append("")
if isinstance(test.declaration, BootstrapMeasurementError):
variables = ", ".join(test.declaration.variables)
lines.append(
f"{test.name}: variables={variables}; resampling={test.resampling}; "
f"original={_format_number(test.declared_effect)}; "
f"successful={len(test.records)}; failed={test.n_failed}"
)
else:
lines.append(
f"{test.name}: process={test.declaration.name}; family={test.family}; "
f"effect={_format_number(test.declared_effect)}; "
f"successful={len(test.records)}; failed={test.n_failed}"
)
lines.append(f"rule: {test.detail}")
if not test.passed and test.declaration is None:
lines.append("")
lines.append(f"{test.name}: {test.detail}")
if self.passed:
lines.append("")
lines.append("VERDICT: all refutation tests behaved as they should.")
return "\n".join(lines)
def __repr__(self) -> str: # pragma: no cover - debugging aid
return self.summary()
_KNOWN_TESTS = (
"placebo",
"random_common_cause",
"subset",
"negative_control_outcome",
"dummy_outcome",
"simulated_outcome",
"bootstrap_measurement_error",
)
_GENERATED_TESTS = ("dummy_outcome", "simulated_outcome")
#: The refutations whose refit sees a different row set from the fit's own. ``subset`` drops
#: rows, so the count changes; ``bootstrap_measurement_error`` resamples with replacement
#: through :meth:`~cleverly.inference.bootstrap._BootstrapDesign.sample`, so the count is
#: unchanged but a different unit sits at each position -- the same reason the targeted
#: bootstrap refuses a supplied plan outright. Every other test here replaces or perturbs a
#: column and keeps each row in its position, which is what makes a supplied
#: :class:`~cleverly.SplitPlan` still apply.
_ROW_SET_TESTS = ("subset", "bootstrap_measurement_error")
_CHILD_SEED_TAGS = {"dummy_outcome": 1, "simulated_outcome": 2}
_ADDITIVE_MEAN_CONTRASTS = {"ate", "att", "atc"}
def _permuted_treatment(data: Any, rng: np.random.Generator) -> np.ndarray:
"""The placebo treatment: the recorded codes permuted among the recorded rows.
A vector treatment with a continuous component permutes its component rows. A row
with an unrecorded treatment keeps its ``NaN`` code, so the placebo keeps the treatment
observation pattern and permutes only what was recorded. Without a missing treatment
this is a permutation of the whole column, drawn exactly as before.
"""
if data.treatment_kind == "vector":
# The rows of the component block move together, so the placebo keeps the joint
# law of the components and breaks only their link to the covariates and outcome.
return np.asarray(rng.permutation(data.treatment_components, axis=0), dtype=float)
if not data.has_missing_treatment:
return np.asarray(rng.permutation(data.treatment), dtype=float)
recorded = data.treatment_recorded
permuted = np.asarray(data.treatment, dtype=float).copy()
permuted[recorded] = rng.permutation(permuted[recorded])
return permuted
def _validated_tests(
tests: Sequence[str], n_replicates: int | None, negative_control_outcome: Any = None
) -> tuple[str, ...]:
"""Refuse malformed arguments of :func:`refute` before any refit.
A malformed argument is reported before any refusal, so a caller who misspelled a test
name hears that before a refusal that the corrected name might not meet. :func:`refute`
and the ``refute`` method of :class:`~cleverly.assessment.DiagnosticsFacade` each call
this first, and :func:`refute_refusal` returns ``None`` where it raises.
Parameters
----------
tests : sequence of str
The requested test names.
n_replicates : int or None
Replicates per randomized test.
negative_control_outcome : Any
The outcome required by the negative-control test.
Returns
-------
tuple of str
The requested test names.
Raises
------
ValueError
If a test name is unknown, ``n_replicates`` is not a positive integer, or the
negative-control test has no outcome.
"""
requested = tuple(tests)
unknown = [name for name in requested if name not in _KNOWN_TESTS]
if unknown:
raise ValueError(
f"unknown refutation test {unknown[0]!r}; choose from {list(_KNOWN_TESTS)}"
)
if n_replicates is not None and (
isinstance(n_replicates, bool)
or not isinstance(n_replicates, (int, np.integer))
or n_replicates < 1
):
raise ValueError("n_replicates must be positive and an integer")
if "negative_control_outcome" in requested and negative_control_outcome is None:
raise ValueError(
"the negative_control_outcome test needs an outcome array that treatment "
"cannot affect; pass negative_control_outcome=<array>"
)
return requested
@dataclass(frozen=True)
class _RefuteRequest:
"""The arguments of one :func:`refute` request that a rule of the table reads.
Parameters
----------
estimand : str
The alias the tests would refit.
tests : tuple of str
The requested test names, already validated.
n_replicates : int or None
Replicates per randomized test.
dummy_outcome : GaussianIndependentOutcome or None
The independent process declaration, or ``None`` for the default.
simulated_outcome : GaussianAdjustmentOutcome or None
The adjustment-dependent process declaration, or ``None`` for the default.
outcome_rule : EmpiricalInclusionRule
The rule applied to generated-outcome refits.
bootstrap_measurement_error : BootstrapMeasurementError or None
The measurement-error declaration.
measurement_error_rule : EmpiricalInclusionRule
The rule applied to measurement-error refits.
"""
estimand: str
tests: tuple[str, ...]
n_replicates: int | None = None
dummy_outcome: Any = None
simulated_outcome: Any = None
outcome_rule: Any = field(default_factory=EmpiricalInclusionRule)
bootstrap_measurement_error: Any = None
measurement_error_rule: Any = field(default_factory=EmpiricalInclusionRule)
@property
def generated(self) -> tuple[str, ...]:
"""The requested generated-outcome tests, in request order."""
return tuple(name for name in self.tests if name in _GENERATED_TESTS)
@property
def draws(self) -> int:
"""The draw budget of a generated-outcome or measurement-error test."""
return DEFAULT_OUTCOME_REPLICATES if self.n_replicates is None else self.n_replicates
def process(self, name: str) -> Any:
"""The outcome process a generated-outcome test draws from.
Parameters
----------
name : str
``"dummy_outcome"`` or ``"simulated_outcome"``.
Returns
-------
GaussianIndependentOutcome or GaussianAdjustmentOutcome
The declared process, or the default one when the request declares none.
"""
if name == "dummy_outcome":
return (
GaussianIndependentOutcome() if self.dummy_outcome is None else self.dummy_outcome
)
return (
GaussianAdjustmentOutcome()
if self.simulated_outcome is None
else self.simulated_outcome
)
#: What a caller can run when a generated-outcome test has no derivation for the fit.
_STABILITY_ALTERNATIVE = (
"The stability tests run on this fit: refute(result, estimand=..., "
"tests=('random_common_cause', 'subset'))."
)
def _invalid(text: str) -> Refusal:
"""The refusal of a request whose declaration the test cannot run on."""
return Refusal(RefusalReason.INVALID_DECLARATION, text)
def _refuse_missing_estimator(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse a result that no longer holds the estimator every test refits."""
if result.estimator is not None:
return None
return _invalid("refute needs the fitted estimator that produced the result")
def _refuse_learned_rule(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse a learned-rule fit, whose refit relearns the rules and so its target."""
from ..interventions.learned import LEARNED_RULE_REFIT_REFUSAL
# Read defensively: a stub result without a configuration meets the later rules.
if getattr(getattr(result, "config", None), "parameter_axis", None) != "learned_rule":
return None
return LEARNED_RULE_REFIT_REFUSAL
def _refuse_unreported_estimand(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse an estimand the fit did not report, which no refit can compare against."""
if request.estimand in result.estimates:
return None
return _invalid(f"estimand {request.estimand!r} was not requested in this fit")
def _target_declaration(result: TMLEResult, estimand: str) -> Any:
"""Read the registered target through a key or the reported alias convention."""
from ..targets import TARGETS
from ..targets.base import parameter_stem
keys = getattr(result, "parameter_keys", {})
key = keys.get(estimand) if keys else None
target = getattr(key, "estimand", parameter_stem(estimand))
return TARGETS.get(target)
def _no_effect_null(result: TMLEResult, estimand: str) -> float | None:
"""Return the reported-scale no-effect value of a registered parameter, if defined."""
declaration = _target_declaration(result, estimand)
if declaration is None:
return None
if declaration.parameter_axis == "msm":
return None
if declaration.scale == "ratio":
return 1.0
if declaration.scale == "difference" or declaration.name == "paf":
return 0.0
return None
def _no_effect_difference(value: float, null: float) -> float:
"""Return the signed difference on the estimate's inference scale."""
if null == 1.0:
return float(np.log(value))
return value - null
def _stability_value(value: float, is_ratio: bool) -> float:
"""Return a refit value on the scale of its plug-in standard error."""
return float(np.log(value)) if is_ratio else value
def _standardized_distance(distance: float, std_error: float) -> float:
"""Express a null distance in standard-error units, including zero spread."""
if std_error == 0.0:
return 0.0 if distance == 0.0 else float("inf")
return distance / std_error
def _refuse_no_effect_null(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse null-effect tests on parameters with no fixed reported-scale null."""
chosen = [name for name in request.tests if name in ("placebo", "negative_control_outcome")]
if not chosen or _no_effect_null(result, request.estimand) is not None:
return None
keys = getattr(result, "parameter_keys", {})
key = keys.get(request.estimand) if keys else None
from ..targets.base import parameter_stem
target = getattr(key, "estimand", parameter_stem(request.estimand))
label = "NaturalCourseMean" if target == "ey_obs" else repr(request.estimand)
prefix = "the placebo refutation and " if target == "ey_obs" and "placebo" in chosen else ""
return _invalid(
f"{prefix}refutation test(s) {chosen} require a fixed no-effect value on the reported "
f"scale, but {label} has none. An outcome level or policy mean need not approach "
"zero, and an arbitrary MSM coefficient has no universal null. Drop these "
"tests from tests=; stability tests keep their usual interpretation."
)
def _refuse_row_set_under_plan(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse a test that refits on another row set when the fit was given a split plan.
A fit that was given a :class:`~cleverly.SplitPlan` cannot refit on a different row
set: the plan holds one fold label per row of the original data, and there is no rule
here for which labels a resampled row inherits. Left to the refit, ``subset`` gives the
caller a row-count ``DataError`` naming a "split plan" they did not mention in this
call, after paying for the tests that ran first, and ``bootstrap_measurement_error``
gives no error at all: its draw keeps the row count, so the labels land on resampled
units and the refutation silently reports a split the fit never ran.
``getattr`` because ``TMLEResult.estimator`` is annotated ``Any`` and :func:`refute`
reads it as a refit *seam*: anything with the ``refit`` signature serves, and
``_validate_generated_eligibility`` reads ``getattr(estimator, "family", "auto")`` off
the same object for the same reason. The stand-ins in
``tests/unit/test_generated_outcome_refutation.py`` and
``tests/unit/test_bootstrap_measurement_error.py`` hold no plan.
"""
supplied_plan = getattr(result.estimator, "split_plan", None)
if supplied_plan is None:
return None
resampled = [name for name in request.tests if name in _ROW_SET_TESTS]
if not resampled:
return None
return _invalid(
f"refutation test(s) {resampled} refit on rows this fit did not run, and "
f"this fit declared split_plan={supplied_plan!r}. A supplied plan labels "
"the rows it was realised on, by position, so it cannot label a refit that "
"drops rows or draws them with replacement, and dropping the plan for those "
"refits would refute a different split from the one the fit ran. Drop "
f"{resampled} from tests=, or refit the result without split_plan= to refute "
"a fit whose folds are drawn from the data each time"
)
def _refuse_generated_rule(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse a generated-outcome test under a rule other than the registered one."""
if not request.generated or type(request.outcome_rule) is EmpiricalInclusionRule:
return None
return _invalid(
"generated-outcome refutations require the exact registered "
"EmpiricalInclusionRule declaration"
)
def _refuse_generated_eligibility(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse a generated-outcome test whose fit, estimand or process has no derivation.
The sentence is the one ``_validate_generated_eligibility`` raises. This rule calls
that check rather than copying it, and :func:`refute` calls it again to read the
contrast direction it returns.
"""
for name in request.generated:
try:
_validate_generated_eligibility(result, request.estimand, request.process(name), name)
except CapabilityError as error:
return error.refusal
return None
def _refuse_measurement_declaration(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse ``bootstrap_measurement_error`` without the registered declaration."""
if "bootstrap_measurement_error" not in request.tests:
return None
if type(request.bootstrap_measurement_error) is BootstrapMeasurementError:
return None
return _invalid(
"bootstrap_measurement_error requires the exact registered "
"BootstrapMeasurementError declaration"
)
def _refuse_measurement_rule(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse ``bootstrap_measurement_error`` under a rule other than the registered one."""
if "bootstrap_measurement_error" not in request.tests:
return None
if type(request.measurement_error_rule) is EmpiricalInclusionRule:
return None
return _invalid(
"bootstrap_measurement_error requires the exact registered "
"EmpiricalInclusionRule declaration"
)
def _refuse_measurement_eligibility(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse a measurement-error declaration that this fit's data cannot carry.
The sentence is the one ``_validate_measurement_error_eligibility`` raises, and this
rule calls that check rather than copying it.
"""
if "bootstrap_measurement_error" not in request.tests:
return None
try:
_validate_measurement_error_eligibility(
result, request.bootstrap_measurement_error, request.draws
)
except CapabilityError as error:
return error.refusal
return None
def _refuse_measurement_budget(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse a measurement-error draw budget below its rule's floor."""
if "bootstrap_measurement_error" not in request.tests:
return None
floor = request.measurement_error_rule.minimum_draws
if request.draws >= floor:
return None
return _invalid(
f"bootstrap_measurement_error was asked for {request.draws} draw(s) "
f"under a rule that requires {floor}; raise "
"n_replicates or declare a rule with a smaller minimum_draws"
)
def _refuse_generated_budget(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""Refuse a generated-outcome draw budget below its rule's floor.
A budget below the rule's own floor can only end in "too few draws", and the caller
can see that before any refit is paid for.
"""
floor = request.outcome_rule.minimum_draws
if not request.generated or request.draws >= floor:
return None
return _invalid(
f"a generated-outcome refutation was asked for {request.draws} draw(s) under a rule "
f"that requires {floor}; raise n_replicates to at least "
f"{floor}, or declare a rule with a smaller minimum_draws"
)
#: Every refusal of a :func:`refute` request that its arguments and the result decide
#: before any refit, in the one order the call and the ``refute`` capability row both use.
#: Each rule assumes its predecessors returned ``None``. The names are the introspection
#: contract, so a test pins the order without respelling a message.
_REQUEST_RULES: tuple[tuple[str, Callable[[TMLEResult, _RefuteRequest], Refusal | None]], ...] = (
("estimator", _refuse_missing_estimator),
("learned_rule", _refuse_learned_rule),
("estimand", _refuse_unreported_estimand),
("no_effect_null", _refuse_no_effect_null),
("row_set_under_plan", _refuse_row_set_under_plan),
("generated_rule", _refuse_generated_rule),
("generated_eligibility", _refuse_generated_eligibility),
("measurement_declaration", _refuse_measurement_declaration),
("measurement_rule", _refuse_measurement_rule),
("measurement_eligibility", _refuse_measurement_eligibility),
("measurement_budget", _refuse_measurement_budget),
("generated_budget", _refuse_generated_budget),
)
def _first_refusal(result: TMLEResult, request: _RefuteRequest) -> Refusal | None:
"""The record of the first rule of :data:`_REQUEST_RULES` that refuses ``request``."""
for _name, rule in _REQUEST_RULES:
refusal = rule(result, request)
if refusal is not None:
return refusal
return None
def refute_refusal(
result: TMLEResult,
*,
estimand: str,
tests: Sequence[str],
n_replicates: int | None = None,
negative_control_outcome: Any = None,
dummy_outcome: GaussianIndependentOutcome | None = None,
simulated_outcome: GaussianAdjustmentOutcome | None = None,
outcome_rule: EmpiricalInclusionRule = EmpiricalInclusionRule(),
bootstrap_measurement_error: BootstrapMeasurementError | None = None,
measurement_error_rule: EmpiricalInclusionRule = EmpiricalInclusionRule(),
) -> Refusal | None:
"""Return the refusal a :func:`refute` request meets before any refit.
:func:`refute` raises the record this function returns, and the ``refute`` capability
row of :class:`~cleverly.assessment.DiagnosticsFacade` quotes it for the same request.
The row defers on ``tests`` when the default tests are refused and a smaller set runs,
so a combined report names the argument instead of calling a refutation the call
refuses. The keyword arguments are those of :func:`refute` that decide a refusal.
A malformed argument is not a refusal. An unknown test name, an invalid
``n_replicates``, or a missing negative-control outcome raises ``ValueError`` from
the call before any refusal. This function returns ``None`` for those requests.
Parameters
----------
result : TMLEResult
A fitted point-treatment result.
estimand : str
The alias the tests would refit.
tests : sequence of str
The requested test names.
n_replicates : int or None
Replicates per randomized test.
negative_control_outcome : Any
The outcome required by the negative-control test.
dummy_outcome : GaussianIndependentOutcome or None
Independent Gaussian process declaration.
simulated_outcome : GaussianAdjustmentOutcome or None
Adjustment-dependent additive process declaration.
outcome_rule : EmpiricalInclusionRule
Rule applied to generated-outcome refits.
bootstrap_measurement_error : BootstrapMeasurementError or None
Measurement-error declaration.
measurement_error_rule : EmpiricalInclusionRule
Rule applied to measurement-error refits.
Returns
-------
Refusal or None
The exact refusal record, or ``None`` when no rule refuses the request or an
argument is malformed.
"""
try:
requested = _validated_tests(tests, n_replicates, negative_control_outcome)
except ValueError:
return None
request = _RefuteRequest(
estimand,
requested,
n_replicates,
dummy_outcome,
simulated_outcome,
outcome_rule,
bootstrap_measurement_error,
measurement_error_rule,
)
return _first_refusal(result, request)
def _generated_child_seeds(root_seed: int, name: str, count: int) -> tuple[int, ...]:
sequence = np.random.SeedSequence([root_seed, _CHILD_SEED_TAGS[name]])
return tuple(int(child.generate_state(1)[0]) for child in sequence.spawn(count))
def _classification_only_learner(learner: Any) -> Any | None:
"""Return a saved learner only when it has no regression-capable route."""
if learner is None:
return None
from sklearn.base import is_classifier
if getattr(learner, "task", None) == "classification" or is_classifier(learner):
return learner
library = getattr(learner, "library", None)
if library is None:
return None
incompatible_candidates = []
for item in library:
candidate = item[1] if isinstance(item, tuple) and len(item) == 2 else item
incompatible_candidates.append(_classification_only_learner(candidate) is not None)
return learner if incompatible_candidates and all(incompatible_candidates) else None
def _validate_generated_process(process: Any, name: str) -> None:
expected = GaussianIndependentOutcome if name == "dummy_outcome" else GaussianAdjustmentOutcome
if type(process) is not expected:
raise CapabilityError(
f"{name} has no registered effect derivation and family validation for "
f"{type(process).__name__}; use the exact registered {expected.__name__} declaration",
reason=RefusalReason.INVALID_DECLARATION,
)
if process.family != "gaussian":
raise CapabilityError(
f"{name} has family={process.family!r}; only the Gaussian process family has "
"an implemented effect derivation",
reason=RefusalReason.INVALID_DECLARATION,
)
def _resolve_arm_code(data: Any, endpoint: Any, role: str, name: str) -> float:
"""Resolve one parameter-key endpoint to its internal arm code.
Parameters
----------
data : CausalData
Analysis rows that declare the arm codes and their labels.
endpoint : Any
Label or code recorded on the parameter key.
role : str
Which endpoint this is, named in a refusal.
name : str
Refutation name, named in a refusal.
Returns
-------
float
The internal arm code for ``endpoint``.
Raises
------
CapabilityError
If no declared arm matches the endpoint.
"""
# Same label-or-code match the estimator's reference-arm resolution uses, so a key
# written in the user's own levels and one written in internal codes agree.
for code, label in zip(data.arm_codes, data.treatment_levels, strict=True):
if label == endpoint or code == endpoint:
return float(code)
raise CapabilityError(
f"{name} has no effect derivation for {role}={endpoint!r}; it is not an arm of "
f"{data.treatment_name} whose levels are {list(data.treatment_levels)}",
reason=RefusalReason.INVALID_DECLARATION,
)
def _validate_generated_eligibility(
result: Any, estimand: str, process: Any, name: str
) -> tuple[float, float]:
"""Refuse an ineligible generated-outcome request and resolve its contrast direction.
Parameters
----------
result : TMLEResult
The fitted result the refutation would refit.
estimand : str
Alias the refutation runs for.
process : object
Registered outcome-process declaration.
name : str
Refutation name, named in every refusal.
Returns
-------
tuple of float, float
Value and reference arm codes the parameter key contrasts, in that order.
Raises
------
CapabilityError
If the fit, the estimand, or the process has no registered effect derivation.
"""
_validate_generated_process(process, name)
from ..study import (
BackdoorMeanContrast,
ExplicitAdjustmentProvider,
ParameterKey,
_matches_registered_point_identification,
)
from ..targets import TARGETS
identified = getattr(result, "identified_effect", None)
if identified is None:
raise CapabilityError(
f"{name} needs identification metadata for a backdoor additive mean contrast; "
"this result records none",
reason=RefusalReason.INVALID_DECLARATION,
)
functional = identified.functional
if type(functional) is not BackdoorMeanContrast:
raise CapabilityError(
f"{name} has no outcome process or effect derivation for {type(functional).__name__}.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative="Use a backdoor-identified effect.",
roadmap="X67",
)
provider = getattr(identified, "provider", None)
if type(provider) is not ExplicitAdjustmentProvider:
raise CapabilityError(
f"{name} needs registered backdoor provider provenance; "
f"{type(provider).__name__} does not explicitly provide it",
reason=RefusalReason.INVALID_DECLARATION,
)
data = result.data
if data.family != process.family:
raise CapabilityError(
f"{name} cannot replace an original outcome family={data.family!r} with the "
f"registered {process.family!r} process; the first catalog requires matching "
"Gaussian outcome configurations.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=_STABILITY_ALTERNATIVE,
roadmap="X67",
)
estimator = result.estimator
configured_family = getattr(estimator, "family", "auto")
if configured_family not in {"auto", process.family}:
raise CapabilityError(
f"{name} cannot refit an estimator configured with family={configured_family!r} "
f"under the declared {process.family!r} process",
reason=RefusalReason.INVALID_DECLARATION,
)
if getattr(estimator, "q_bounds", None) is not None:
raise CapabilityError(
f"{name} replaces the outcome with a draw from the registered Gaussian "
"process, whose support is the whole real line, and this fit declares "
f"q_bounds={estimator.q_bounds!r}. The refit would scale the generated "
"outcome into a range it is not confined to and report the refutation under "
"bounds the process does not satisfy. Leave q_bounds=None on the fit you "
"refute (Targeting(q_bounds=None)). A cross-fitted continuous fit has to "
"declare them, so refute a continuous outcome on a fit with cross_fit=False "
"(CrossFitting(enabled=False))",
reason=RefusalReason.INVALID_DECLARATION,
)
outcome_learner = getattr(estimator, "outcome_learner", None)
incompatible_learner = _classification_only_learner(outcome_learner)
if incompatible_learner is not None:
raise CapabilityError(
f"{name} needs a regression-capable saved outcome learner for its Gaussian "
f"process; {type(incompatible_learner).__name__} is configured for classification",
reason=RefusalReason.INVALID_DECLARATION,
)
if getattr(functional, "longitudinal", False):
raise CapabilityError(
f"{name} has no longitudinal outcome process or sequential effect derivation.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=NO_ALTERNATIVE,
roadmap="X67",
)
if not getattr(data, "is_binary_treatment", False):
raise CapabilityError(
f"{name} has no effect derivation for a non-binary treatment; the first "
"process catalog covers code one against code zero.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=_STABILITY_ALTERNATIVE,
roadmap="X67",
)
if getattr(data, "has_intermediate", False) or functional.intermediate is not None:
raise CapabilityError(
f"{name} has no controlled-direct-effect outcome process or effect derivation.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=_STABILITY_ALTERNATIVE,
roadmap="X67",
)
if getattr(data, "has_missing_outcome", False):
raise CapabilityError(
f"{name} has no missing-outcome process or observation-model effect derivation.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=_STABILITY_ALTERNATIVE,
roadmap="X67",
)
if functional.msm is not None or functional.axis == "msm":
raise CapabilityError(
f"{name} has no MSM outcome process or coefficient derivation.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=_STABILITY_ALTERNATIVE,
roadmap="X67",
)
if functional.axis != "arm" or functional.interventions:
raise CapabilityError(
f"{name} has no intervention-indexed outcome process or effect derivation.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=_STABILITY_ALTERNATIVE,
roadmap="X67",
)
key = result.parameter_keys.get(estimand)
if type(key) is not ParameterKey:
raise CapabilityError(
f"{name} needs a structured parameter key for {estimand!r}; choose one of "
f"{list(result.parameter_keys)}",
reason=RefusalReason.INVALID_DECLARATION,
)
if key.estimand not in _ADDITIVE_MEAN_CONTRASTS:
kind = "ratio" if key.estimand in {"rr", "or"} else key.estimand
raise CapabilityError(
f"{name} has no {kind} effect derivation.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative="Choose an additive mean contrast (ate, att, or atc).",
roadmap="X67",
)
declared_estimand = getattr(getattr(identified, "estimand", None), "name", None)
if functional.target != key.estimand or declared_estimand != key.estimand:
raise CapabilityError(
f"{name} found inconsistent registered target provenance: functional target "
f"{functional.target!r}, identified estimand {declared_estimand!r}, and parameter "
f"key {key.estimand!r} must agree",
reason=RefusalReason.INVALID_DECLARATION,
)
registered = TARGETS.get(functional.target)
declared_identification = getattr(identified, "identification", None)
if (
registered is None
or declared_identification is None
or not _matches_registered_point_identification(
declared_identification,
registered.identification,
identified,
data,
functional.axis,
)
):
raise CapabilityError(
f"{name} needs the registered identification artifact for target {functional.target!r}",
reason=RefusalReason.INVALID_DECLARATION,
)
# The key names the contrast direction, and the process declares its effect for code
# one against code zero. A fit reported against a non-default reference therefore
# declares the negated effect, and reading the key is the only way to see it.
reference_endpoint = key.reference
reference_code = (
float(data.arm_codes[0])
if reference_endpoint is None
else _resolve_arm_code(data, reference_endpoint, "the parameter key reference", name)
)
value_endpoint = key.value
if value_endpoint is None:
remaining = [float(code) for code in data.arm_codes if float(code) != reference_code]
if len(remaining) != 1:
raise CapabilityError(
f"{name} has no effect derivation for a parameter key with no value arm; "
f"{data.treatment_name} does not leave exactly one other arm",
reason=RefusalReason.INVALID_DECLARATION,
)
value_code = remaining[0]
else:
value_code = _resolve_arm_code(data, value_endpoint, "the parameter key value", name)
if value_code == reference_code:
raise CapabilityError(
f"{name} has no effect derivation for a parameter key whose value and reference "
f"are the same arm ({value_code:g}); the catalog covers a two-arm contrast",
reason=RefusalReason.INVALID_DECLARATION,
)
return value_code, reference_code
def _run_empirical_refits(
result: Any,
*,
estimand: str,
draws: Sequence[tuple[int, int]],
replacement: Callable[[int, int], Any],
expected_family: str,
retain_preparation_failures: bool = False,
) -> tuple[tuple[EmpiricalRefitRecord, ...], tuple[ReplicationFailure, ...]]:
"""Run, validate, and retain empirical refits through one failure contract.
Parameters
----------
result : TMLEResult
Fitted result that supplies the estimator and its intermediate value.
estimand : str
Structured result alias read off each refit.
draws : sequence of tuple of int
One ``(replicate, seed)`` pair per requested draw.
replacement : callable
Builds the data for one draw from its replicate index and seed.
expected_family : str
Outcome family every successful refit must report.
retain_preparation_failures : bool
Whether a failure inside ``replacement`` becomes a retained failure. ``False``
raises it, because a deterministic preparation fails the same way on every draw
and the caller must see that declaration or data error at once. ``True`` retains
it, because a bootstrap draw fails per sample: a resample can hold too few rows,
and ``run_bootstrap`` drops exactly that draw. Keep the two paths apart. A refit
failure is retained either way.
Returns
-------
tuple of tuple
Successful records, then the retained failures.
"""
records: list[EmpiricalRefitRecord] = []
failures: list[ReplicationFailure] = []
for replicate, child_seed in draws:
try:
prepared = replacement(replicate, child_seed)
except Exception as error:
if not retain_preparation_failures:
raise
# A bootstrap sample this draw cannot build, or a perturbation the sample
# refuses, is one failed draw among many. Raising would discard every draw
# that already succeeded, so record it with the fields a refit failure carries.
failures.append(
ReplicationFailure(
replicate=replicate,
seed=child_seed,
error_type=type(error).__name__,
message=str(error),
)
)
continue
try:
refitted = result.estimator.refit(
prepared,
intermediate_value=result.intermediate_value,
random_state=child_seed,
)
except Exception as error: # an estimator failure belongs in the report
failures.append(
ReplicationFailure(
replicate=replicate,
seed=child_seed,
error_type=type(error).__name__,
message=str(error),
)
)
continue
fitted_family = str(refitted.data.family)
if fitted_family != expected_family:
raise RuntimeError(
f"empirical refit reported authoritative family={fitted_family!r}; expected "
f"{expected_family!r}"
)
estimate = refitted[estimand]
if not np.isfinite(estimate.psi) or not np.isfinite(estimate.plugin_std_error):
failures.append(
ReplicationFailure(
replicate=replicate,
seed=child_seed,
error_type="ValueError",
message="refitted estimate or standard error is non-finite",
)
)
continue
records.append(
EmpiricalRefitRecord(
replicate=replicate,
seed=child_seed,
estimand=estimand,
estimate=float(estimate.psi),
std_error=float(estimate.plugin_std_error),
family=fitted_family,
)
)
return tuple(records), tuple(failures)
def _measurement_encodings(data: Any) -> dict[str, Any]:
return {encoding.column: encoding for encoding in data.encodings}
def _validate_numeric_measurement_column(data: Any, name: str) -> None:
"""Refuse a selected column that relative Gaussian noise cannot perturb coherently.
Parameters
----------
data : CausalData
Original data the fit ran on.
name : str
Selected adjustment variable that carries no ``CategoricalEncoding``.
"""
column = data.covariate_names.index(name)
values = np.asarray(data.covariates[:, column], dtype=float)
# A column with no encoding is treated as numeric, and ``CausalData.from_arrays``
# records no encoding at all. An undeclared 0/1 column therefore reaches the numeric
# path, where continuous noise would make it real-valued, and a hand-built dummy block
# would have each indicator perturbed on its own. Refuse both here, before any refit.
if bool(np.all(np.isin(values, (0.0, 1.0)))):
raise CapabilityError(
f"bootstrap_measurement_error cannot perturb {name!r}: its values are all 0 or 1 "
"but it carries no CategoricalEncoding, so the operation cannot tell an "
"undeclared binary variable from one indicator of a dummy block, and relative "
"Gaussian noise would make it real-valued. Supply the variable as a boolean or "
"categorical dataframe column to CausalData.from_frame, so it carries the "
"logical metadata this operation perturbs by",
reason=RefusalReason.INVALID_DECLARATION,
)
if float(np.std(values, ddof=0)) == 0.0:
raise CapabilityError(
f"bootstrap_measurement_error cannot perturb {name!r}: the variable is constant, "
"so its relative noise scale is exactly zero and every draw would return the "
"original column while the report still named the variable as perturbed. Remove "
"it from the declaration, or declare an absolute noise law for it",
reason=RefusalReason.INVALID_DECLARATION,
)
def _validate_measurement_error_eligibility(
result: Any, declaration: BootstrapMeasurementError, n_replicates: int
) -> Any:
"""Validate the complete bootstrap measurement-error request before refitting."""
from ..data import CausalData
data = getattr(result, "data", None)
if type(data) is not CausalData:
raise CapabilityError(
"bootstrap_measurement_error supports point-treatment CausalData results; "
f"got result family {type(result).__name__} with data family {type(data).__name__}.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=NO_ALTERNATIVE,
roadmap="X67",
)
if n_replicates < 1:
raise ValueError("n_replicates must be positive")
if declaration.resampling == "cluster" and data.cluster is None:
raise CapabilityError(
"resampling='cluster' requires the data to carry cluster ids",
reason=RefusalReason.INVALID_DECLARATION,
)
selected_strata = set(declaration.variables).intersection(data.strata_names)
if selected_strata:
raise CapabilityError(
"bootstrap_measurement_error cannot perturb selected strata variables while "
f"preserving target metadata: {sorted(selected_strata)}.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=(
"Declare BootstrapMeasurementError(variables=...) with the adjustment "
"variables other than the strata."
),
roadmap="X67",
)
encodings = _measurement_encodings(data)
generated = data.indicator_sources
numeric = set(data.covariate_names).difference(generated)
original = numeric.union(encodings)
unknown = [name for name in declaration.variables if name not in original]
if unknown:
indicator = next((name for name in unknown if name in generated), None)
if indicator is not None:
raise CapabilityError(
f"{indicator!r} is a generated indicator for original categorical variable "
f"{generated[indicator]!r}.",
reason=RefusalReason.API_REDIRECTION,
alternative=f"Select the original variable {generated[indicator]!r}.",
)
raise CapabilityError(
f"unknown adjustment variables {unknown}; choose from {sorted(original)}",
reason=RefusalReason.INVALID_DECLARATION,
)
retained = set(data.covariate_names)
for name in declaration.variables:
encoding = encodings.get(name)
if encoding is None:
_validate_numeric_measurement_column(data, name)
continue
missing = [indicator for indicator in encoding.generated if indicator not in retained]
if missing:
raise CapabilityError(
"bootstrap_measurement_error cannot perturb categorical variable "
f"{name!r} because duplicate-column removal left its encoded block "
f"incomplete; missing indicators {missing}",
reason=RefusalReason.INVALID_DECLARATION,
)
return data
def _categorical_codes(data: Any, encoding: Any) -> np.ndarray:
columns = [data.covariate_names.index(name) for name in encoding.generated]
block = np.asarray(data.covariates[:, columns], dtype=float)
valid = np.all(np.isin(block, (0.0, 1.0)), axis=1) & (np.sum(block, axis=1) <= 1.0)
if not np.all(valid):
raise RuntimeError(
f"encoded block for {encoding.column!r} is not a valid drop-first indicator block"
)
active = np.argmax(block, axis=1) + 1
return np.where(np.sum(block, axis=1) == 0.0, 0, active).astype(np.int64)
def _perturb_measurement_error(
data: Any, declaration: BootstrapMeasurementError, *, seed: int
) -> Any:
"""Perturb sampled adjustment variables and rebuild complete categorical blocks."""
rng = np.random.default_rng(seed)
values = np.array(data.covariates, dtype=float, copy=True)
encodings = _measurement_encodings(data)
for name in declaration.variables:
encoding = encodings.get(name)
if encoding is None:
column = data.covariate_names.index(name)
values[:, column] += declaration.numeric_noise.draw(rng, values[:, column])
continue
columns = [data.covariate_names.index(item) for item in encoding.generated]
codes = _categorical_codes(data, encoding)
change = rng.random(data.n) < declaration.categorical_change_probability
if np.any(change):
alternatives = rng.integers(0, len(encoding.levels) - 1, size=int(np.sum(change)))
current = codes[change]
codes[change] = alternatives + (alternatives >= current)
rebuilt = np.column_stack(
[np.asarray(codes == code, dtype=float) for code in range(1, len(encoding.levels))]
)
values[:, columns] = rebuilt
return data.with_covariates(values, name="measurement-error covariates")
def _bootstrap_measurement_error_test(
result: Any,
*,
estimand: str,
original: float,
declaration: BootstrapMeasurementError,
rule: EmpiricalInclusionRule,
n_replicates: int,
root_seed: int,
) -> RefutationTest:
data = _validate_measurement_error_eligibility(result, declaration, n_replicates)
design = _bootstrap_design(
data,
n_replicates=n_replicates,
resampling=declaration.resampling,
random_state=root_seed,
)
draw_by_replicate = {draw.replicate: draw for draw in design.draws}
records, failures = _run_empirical_refits(
result,
estimand=estimand,
draws=tuple((draw.replicate, draw.seed) for draw in design.draws),
replacement=lambda replicate, seed: _perturb_measurement_error(
design.sample(data, draw_by_replicate[replicate]), declaration, seed=seed
),
expected_family=data.family,
# A cluster or stratum bootstrap draw can hold too few rows to subset, and that is
# a property of the draw rather than of the request. ``run_bootstrap`` drops such a
# draw, so this operation retains it.
retain_preparation_failures=True,
)
values = tuple(item.estimate for item in records)
passed, probability, detail = rule.evaluate(values, original, failures)
families = {item.family for item in records}
family = next(iter(families)) if len(families) == 1 else None
if len(families) > 1:
passed = False
detail = f"successful refits reported inconsistent outcome families {sorted(families)}"
return RefutationTest(
name="bootstrap_measurement_error",
estimand=estimand,
original=original,
values=values,
expectation=f"includes {original:.5g}",
passed=passed,
detail=detail,
standard_errors=tuple(item.std_error for item in records),
declaration=declaration,
rule=rule,
records=records,
failures=failures,
family=family,
empirical_pvalue=probability,
declared_effect=original,
requested_draws=n_replicates,
resampling=design.resampling,
)
def _generated_outcome_test(
result: Any,
*,
name: str,
estimand: str,
original: float,
process: GaussianIndependentOutcome | GaussianAdjustmentOutcome,
rule: EmpiricalInclusionRule,
n_replicates: int,
root_seed: int,
value_code: float,
reference_code: float,
) -> RefutationTest:
data = result.data
seeds = _generated_child_seeds(root_seed, name, n_replicates)
records, failures = _run_empirical_refits(
result,
estimand=estimand,
draws=tuple(enumerate(seeds)),
replacement=lambda _replicate, seed: data.with_outcome(
process.draw(data, seed=seed), family=process.family
),
expected_family=process.family,
)
values = tuple(item.estimate for item in records)
# The additive treatment term makes the arm contrast linear in the code difference, so
# the declared effect is the process effect signed by the key's own direction.
declared_effect = float(process.known_effect * (value_code - reference_code))
passed, probability, detail = rule.evaluate(values, declared_effect, failures)
families = {item.family for item in records}
family = next(iter(families)) if len(families) == 1 else None
if len(families) > 1:
passed = False
detail = f"successful refits reported inconsistent outcome families {sorted(families)}"
return RefutationTest(
name=name,
estimand=estimand,
original=original,
values=values,
expectation=f"includes {declared_effect:.5g}",
passed=passed,
detail=detail,
standard_errors=tuple(item.std_error for item in records),
declaration=process,
rule=rule,
records=records,
failures=failures,
family=family,
empirical_pvalue=probability,
declared_effect=declared_effect,
requested_draws=n_replicates,
)
[docs]
def refute(
result: TMLEResult,
*,
estimand: str = "ate",
tests: Sequence[str] = DEFAULT_TESTS,
n_replicates: int | None = None,
subset_fraction: float = 0.7,
negative_control_outcome: Any = None,
dummy_outcome: GaussianIndependentOutcome | None = None,
simulated_outcome: GaussianAdjustmentOutcome | None = None,
outcome_rule: EmpiricalInclusionRule = EmpiricalInclusionRule(),
bootstrap_measurement_error: BootstrapMeasurementError | None = None,
measurement_error_rule: EmpiricalInclusionRule = EmpiricalInclusionRule(),
random_state: int | None = None,
tolerance: float = 3.0,
) -> RefutationResult:
"""Run refutation tests against a fitted model.
Parameters
----------
result : TMLEResult
A fitted result to refute.
estimand : str
Alias the tests are run for.
tests : sequence of str
Which refutations to run. Parameters without a fixed reported-scale no-effect
value refuse ``placebo`` and ``negative_control_outcome`` before any refit.
Differences use zero; risk and odds ratios use one. ``subset`` and
``random_common_cause`` keep their stability interpretations. A fit
given ``split_plan=`` refuses ``subset`` and ``bootstrap_measurement_error``
before any refit. ``refute_refusal`` holds these refusals, and every refusal of a
generated-outcome or measurement-error request. The ``refute`` capability row
reads it for the same request, and defers on ``tests`` when the default tests meet
one and a smaller set runs. ``random_common_cause`` adds a column of independent
noise. A :class:`~cleverly.CTMLE` with an explicit ``ordering=`` places that column
after the declared ordering, and a :class:`~cleverly.DRTMLE` refits without an
``evaluation=`` companion that lacks it.
n_replicates : int or None
Replicates per randomized test. ``None`` uses five for each established
perturbation, and 100 for each generated-outcome test and for
``bootstrap_measurement_error``. Each of those three tests refuses a budget below
its rule's ``minimum_draws`` before any refit.
subset_fraction : float
Share of rows the subset test refits on.
negative_control_outcome : str or None
An outcome the treatment cannot affect. Required to run that test.
dummy_outcome : GaussianIndependentOutcome or None
Independent Gaussian process declaration. ``None`` uses the default declaration.
simulated_outcome : GaussianAdjustmentOutcome or None
Adjustment-dependent additive process. ``None`` uses the default declaration.
outcome_rule : EmpiricalInclusionRule
Recorded rule applied to generated-outcome refits.
bootstrap_measurement_error : BootstrapMeasurementError or None
Measurement-error declaration. Required when the named test is requested.
measurement_error_rule : EmpiricalInclusionRule
Recorded rule that compares the original estimate with measurement-error refits.
random_state : int or None
Seed for the randomised tests. ``None`` uses the seed the fit was run with, so a
seeded fit gives the same refutation every time. An unseeded fit draws a seed, and
the refits run under it too. Either way the report records it under
``random_state``, and passing that value back repeats the report.
tolerance : float
How many standard errors a null test may deviate before failing. The default of
3 keeps the false-alarm rate low across several tests.
Returns
-------
RefutationResult
One record per test run, with what it expected and what it saw.
Raises
------
ValueError
If a test name is unknown, ``n_replicates`` is not a positive integer, or the
negative-control test has no outcome. These checks precede every refusal and refit.
CapabilityError
If ``refute_refusal`` refuses the request. Each such refusal precedes every refit.
"""
# The malformed arguments first, and the refusals of a well-formed request after them.
# Every refusal is decided before any refit, so a mixed call does not pay for placebo
# fits and only then discover that its outcome process has no effect derivation. The
# ``refute`` capability row reads the same table, so a row that reads available is a
# request this line admits.
request = _RefuteRequest(
estimand,
_validated_tests(tests, n_replicates, negative_control_outcome),
n_replicates,
dummy_outcome,
simulated_outcome,
outcome_rule,
bootstrap_measurement_error,
measurement_error_rule,
)
refusal = _first_refusal(result, request)
if refusal is not None:
raise CapabilityError.from_refusal(refusal)
requested = request.tests
estimator = result.estimator
processes = {name: request.process(name) for name in _GENERATED_TESTS}
# The table admitted each generated test, so this reads the contrast direction and
# refuses nothing.
contrasts = {
name: _validate_generated_eligibility(result, estimand, processes[name], name)
for name in request.generated
}
original = result[estimand].psi
no_effect_null = _no_effect_null(result, estimand)
estimate_scale = getattr(result[estimand], "scale", None)
declaration = _target_declaration(result, estimand)
is_ratio = (estimate_scale or getattr(declaration, "scale", None)) == "ratio"
# A tolerance scale for the refutation comparisons below, not a coverage claim, so
# the plug-in accessor answers on a selector-path collaborative fit too.
std_error = result[estimand].plugin_std_error
# Every refit reruns this estimator, so each test's standard errors carry the fit's
# own inference status. ``RefutationTest.to_frame`` names its column from it, and a
# detail below names the unit it counts in from it.
status = result[estimand].inference
unit = spread_name("standard errors", status)
# The package convention for a stochastic operation on a fitted object: an explicit seed
# wins, ``None`` inherits the fit's own. ``is None`` rather than truthiness, because
# ``random_state=0`` is falsy and still has to win. Same form as ``ctmle.py:846``.
#
# An unseeded fit draws one seed here rather than handing ``None`` to the generator, so
# the report can name it. The seed then goes to the refits as well as to the draws
# below. A perturbation is only half of a refutation: ``refit`` re-learns the
# nuisances, and an estimator with no ``random_state`` redraws its folds every time, so
# seeding the draws alone would leave the report unrepeatable.
seed = resolve_assessment_seed(result, random_state)
rng = np.random.default_rng(seed)
data = result.data
outcomes: list[RefutationTest] = []
def refit(replacement: Any) -> float:
refitted = estimator.refit(
replacement,
intermediate_value=result.intermediate_value,
random_state=seed,
)
return refitted[estimand].psi
for name in requested:
replicate_count = (
n_replicates
if n_replicates is not None
else (
DEFAULT_OUTCOME_REPLICATES
if name in _GENERATED_TESTS or name == "bootstrap_measurement_error"
else _DEFAULT_PERTURBATION_REPLICATES
)
)
if name == "placebo":
assert no_effect_null is not None
values = tuple(
refit(data.with_treatment(_permuted_treatment(data, rng)))
for _ in range(replicate_count)
)
# Permuting treatment removes the effect, so each replicate is a draw from a
# null distribution whose spread is about the original standard error.
threshold = tolerance * std_error / np.sqrt(max(1, replicate_count))
mean_kind: Literal["arithmetic", "geometric"] = (
"geometric" if is_ratio else "arithmetic"
)
mean = float(np.exp(np.mean(np.log(values)))) if is_ratio else float(np.mean(values))
distance = abs(_no_effect_difference(mean, no_effect_null))
passed = bool(distance <= max(threshold, tolerance * std_error / 2))
outcomes.append(
RefutationTest(
name=name,
estimand=estimand,
original=original,
values=values,
expectation=f"~ {no_effect_null:g}",
passed=passed,
detail=(
f"{mean_kind} mean placebo estimate {mean:+.5g} is "
f"{_standardized_distance(distance, std_error):.1f} {unit} "
f"from {no_effect_null:g}. A large gap can indicate fold leakage "
"or a misdefined estimand"
+ (
". The known treatment mechanism was dropped for each permuted "
"treatment and re-estimated"
if data.known_treatment is not None
else ""
)
),
mean_kind=mean_kind,
)
)
elif name == "random_common_cause":
values = tuple(
refit(data.with_extra_covariate(rng.normal(size=data.n), f"_noise_{index}"))
for index in range(replicate_count)
)
shift = float(np.mean([_stability_value(value, is_ratio) for value in values]))
shift -= _stability_value(original, is_ratio)
passed = bool(abs(shift) <= tolerance * std_error)
outcomes.append(
RefutationTest(
name=name,
estimand=estimand,
original=original,
values=values,
expectation=f"~ {original:.4g}",
passed=passed,
detail=(
f"adding an irrelevant covariate moved the estimate by {shift:+.5g} "
f"{'on the log scale ' if is_ratio else ''}"
f"({_standardized_distance(abs(shift), std_error):.1f} {unit}); "
"a large shift can indicate unstable nuisance models"
),
mean_kind="geometric" if is_ratio else "arithmetic",
)
)
elif name == "subset":
size = max(20, int(subset_fraction * data.n))
values = tuple(
refit(data.subset(rng.choice(data.n, size=size, replace=False)))
for _ in range(replicate_count)
)
expected = std_error * np.sqrt(1.0 / subset_fraction - 1.0 + 1.0)
scaled_values = [_stability_value(value, is_ratio) for value in values]
spread = float(np.std(scaled_values, ddof=1)) if len(values) > 1 else 0.0
passed = bool(spread <= tolerance * expected)
outcomes.append(
RefutationTest(
name=name,
estimand=estimand,
original=original,
values=values,
expectation=(
f"log-scale scatter <~ {tolerance * expected:.4g}"
if is_ratio
else f"scatter <~ {tolerance * expected:.4g}"
),
passed=passed,
detail=(
f"subsample estimates scatter by {spread:.5g} "
f"{'on the log scale ' if is_ratio else ''}against {expected:.5g} "
"expected from sampling alone. A large spread can indicate "
"influential observations"
),
spread_kind="log" if is_ratio else "reported",
)
)
elif name == "negative_control_outcome":
assert no_effect_null is not None
replacement = data.with_outcome(
negative_control_outcome, name="negative_control_outcome"
)
refitted = estimator.refit(
replacement,
intermediate_value=result.intermediate_value,
random_state=seed,
)
value = refitted[estimand].psi
control_se = refitted[estimand].plugin_std_error
distance = abs(_no_effect_difference(value, no_effect_null))
passed = bool(distance <= tolerance * control_se)
outcomes.append(
RefutationTest(
name=name,
estimand=estimand,
original=original,
values=(value,),
expectation=f"~ {no_effect_null:g}",
passed=passed,
detail=(
f"the negative-control outcome shows an estimate of {value:+.5g} "
f"({_standardized_distance(distance, control_se):.1f} {unit} "
f"from the no-effect value {no_effect_null:g}). Under a valid, "
"comparable control design, a large gap can flag residual bias. "
"It can also indicate that the negative-control assumptions fail"
),
)
)
elif name in _GENERATED_TESTS:
outcomes.append(
_generated_outcome_test(
result,
name=name,
estimand=estimand,
original=original,
process=processes[name],
rule=outcome_rule,
n_replicates=replicate_count,
root_seed=seed,
value_code=contrasts[name][0],
reference_code=contrasts[name][1],
)
)
else:
# The preflight above proves the declaration's exact type.
assert bootstrap_measurement_error is not None
outcomes.append(
_bootstrap_measurement_error_test(
result,
estimand=estimand,
original=original,
declaration=bootstrap_measurement_error,
rule=measurement_error_rule,
n_replicates=replicate_count,
root_seed=seed,
)
)
return RefutationResult(
tests=tuple(replace(test, inference=status) for test in outcomes),
estimand=estimand,
backend=result.data.backend,
random_state=seed,
)