r"""Positivity and overlap diagnostics.
Positivity -- every unit having some chance of either treatment given ``W`` -- is the
assumption that fails quietly. Unlike confounding, it leaves a visible fingerprint
in the estimated propensity scores, and unlike model misspecification it is not
fixed by a better learner: if no treated unit resembles a given control unit, no
estimator can say what would have happened to that control unit under treatment.
What to look at, in order of how much it tells you:
**Kish weight concentration.** Kish's equivalent count,
:math:`(\sum_i \omega_i)^2 / \sum_i \omega_i^2` for
:math:`\omega_i = 1/g(W_i)` in the treated arm, describes how concentrated those
nonnegative weights are. A value of 40 out of 500 treated units deserves attention,
but it is not estimator information or the sample size behind the standard error.
**Weight concentration.** The share of total nonnegative diagnostic weight held by the
largest few weights. It excludes the residual, so it is not a share of the realised
estimating equation or influence.
**Truncation load.** How many propensity scores were clipped, and how far. Truncation
trades variance for potential second-order bias. It regularises the finite-sample
procedure, but it does not change the requested estimand. A large clipped fraction
means that the estimate is sensitive to this regularisation in regions with limited
fitted overlap.
Observation weights are folded into :math:`\omega_i` rather than reported separately,
because the diagnostic should describe their joint weight profile with the mechanism
factor. Their separately computed Kish ratios do not in general multiply. For the
weighting profile on its own -- and for the
estimand statement that goes with it -- see
:meth:`~cleverly.data.CausalData.weight_report` and :mod:`cleverly.data.weighting`.
Use :func:`truncation_curve` to see how much the answer actually moves as the
truncation bound changes -- a flat curve is reassuring in a way that no single
diagnostic can be.
"""
from __future__ import annotations
from collections.abc import Callable
from dataclasses import dataclass, field, replace
from typing import TYPE_CHECKING, Any, Literal
import numpy as np
from .._refusal import Refusal, RefusalReason
from .._typing import BoolArray, FloatArray
from ..data.known_mechanism import moved_known_values_refusal
from ..data.validate import MISSING_OUTCOME_DECLARATION
from ..data.weighting import (
REPORTED_DRAW,
SCORE_LOAD_VERDICT_RATIO_FORMAT,
ScoreLoadRow,
effective_sample_size,
score_load_row,
top_weight_share,
validate_score_loads,
)
from ..estimators._nuisance import KnownPropensity
from ..estimators.composite import composite_state
from ..estimators.direct_effect import targeted_rows
from ..estimators.targeting import build_submodel
from ..exceptions import CapabilityError, DataError
from ..inference.influence import median_estimates
from ..targets import TARGETS, parameter_stem
from ..targets.population_intervention import NATURAL_COURSE_SUPPORT_REFUSAL
from ..utils.bounds import g_bounds_for
from ..utils.frames import emit_frame
from ..utils.text import format_draw, format_table
if TYPE_CHECKING: # pragma: no cover - typing only
from ..estimators.base import TMLEResult
from ..fluctuation.submodel import Submodel
__all__ = ["PositivityReport", "positivity_report", "truncation_curve"]
#: Quantiles reported for the propensity distribution in each arm.
_QUANTILES = (0.0, 0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99, 1.0)
#: Thresholds at which the mass of extreme propensity scores is reported.
_THRESHOLDS = (0.01, 0.025, 0.05, 0.1)
class GroupScoreLoadRow(ScoreLoadRow):
"""One equation's descriptive absolute residual-multiplier concentration.
The ten shared concentration keys come from
:class:`~cleverly.data.weighting.ScoreLoadRow`. What a target-group row adds is the
propensity bound that group was held to and how much of the mechanism that bound moved,
neither of which an intervention row has: a regime, shift or tilt is described by the
policy it declares rather than by an arm-level truncation interval.
There is deliberately no draw count here. A group row is held by a
:class:`PositivityReport`, which records :attr:`PositivityReport.n_repeats` once for the
whole fit, and a second copy on every row is a second place for it to disagree with the
first. The intervention rows carry their own because the reports that hold them do not
have such a field.
"""
lower_bound: float
upper_bound: float
clipped_count: float
clipped_fraction: float
#: The derived denominators the clever covariate divides by: ``g_a(W) pi_a(W)`` when the
#: outcome can be missing, ``g_a(W) q_z(a, W)`` for a controlled direct effect, and the
#: three-factor product when a controlled direct effect also has missing outcomes. Each is
#: a *product of separately truncated factors* rather than a fitted or targeted mechanism of
#: its own, which is why they are named once here: the report has to treat them differently
#: from the rows that were held to a single bound.
JOINT_MECHANISM = "P(A=a,Delta=1|W)"
JOINT_INTERMEDIATE = "P(A=a,Z=z|W)"
JOINT_MISSING_INTERMEDIATE = "P(A=a,Delta=1,Z=z|W)"
#: The composite mechanisms of a declared missing treatment, without and with ``delta=``:
#: the probability of the composite indicator ``Delta_A * Delta * 1{A = a}``.
COMPOSITE_TREATMENT = "P(A=a,Delta_A=1|W)"
COMPOSITE_MECHANISM = "P(A=a,Delta_A=1,Delta=1|W)"
#: The treatment observation factor of a declared missing treatment.
TREATMENT_OBSERVATION = "P(Delta_A=1|W)"
#: Which derived row a fit reports, keyed by ``(an observation mechanism was fitted, an
#: intermediate density was, a treatment observation mechanism was)``. The two derived
#: facts below come from this one table, so a further composition is one entry rather than
#: three edits that have to agree. A missing treatment with an intermediate is refused
#: before any fit, so it has no entry.
_COMPOSED_ROWS: dict[tuple[bool, bool, bool], str] = {
(True, False, False): JOINT_MECHANISM,
(False, True, False): JOINT_INTERMEDIATE,
(True, True, False): JOINT_MISSING_INTERMEDIATE,
(False, False, True): COMPOSITE_TREATMENT,
(True, False, True): COMPOSITE_MECHANISM,
}
#: The conditioning sets of a per-arm outcome-adaptive design, without and with baseline
#: strata. Such a fit regresses every factor on ``Qbar(a, W)`` (and ``S``), not on ``W``,
#: so its product is the composite indicator's probability given that design.
_DESIGN_CONDITIONS = ("Qbar(a,W)", "Qbar(a,W),S")
def _on_design(name: str, condition: str) -> str:
"""A derived row's name with the per-arm design in place of ``W`` as the condition."""
return name.replace("|W)", f"|{condition})")
#: The derived rows of a per-arm outcome-adaptive fit, keyed by the row they replace and
#: the design's conditioning set. A controlled direct effect has no per-arm design.
_DESIGN_ROWS: dict[tuple[str, str], str] = {
(name, condition): _on_design(name, condition)
for key, name in _COMPOSED_ROWS.items()
if not key[1]
for condition in _DESIGN_CONDITIONS
}
#: The derived rows of a missing treatment, which also report the share of unit-arm cells
#: below the propensity floor -- the scale R ``drtmle``'s ``tolg`` floors the product at.
_COMPOSITE_ROWS = frozenset(
{COMPOSITE_TREATMENT, COMPOSITE_MECHANISM}
| {
design
for (name, _), design in _DESIGN_ROWS.items()
if name in (COMPOSITE_TREATMENT, COMPOSITE_MECHANISM)
}
)
#: How many ``[nuisance_bound, 1]`` factors stand beside ``g`` in each derived row.
_COMPOSED_NUISANCE_FACTORS = {name: sum(key) for key, name in _COMPOSED_ROWS.items()} | {
design: sum(key)
for key, name in _COMPOSED_ROWS.items()
for (base, _), design in _DESIGN_ROWS.items()
if base == name
}
_COMPOSED_MECHANISMS = frozenset(_COMPOSED_ROWS.values()) | frozenset(_DESIGN_ROWS.values())
#: The target groups whose clever covariate really divides by ``g_a(W) pi_a(W) q_z(a, W)``
#: -- :func:`~cleverly.fluctuation.submodel.mean_submodel`, ``regime_submodel`` and
#: ``msm_submodel`` all build ``1 / (g * pi * pz)`` arm by arm. The rest do not, and a
#: derived row would name a denominator they never form: ``att`` and ``atc`` divide by
#: ``P(A = a)`` in place of ``g_a`` and reweight the reference arm by the propensity odds,
#: ``ipsi`` discards the propensity outright, and ``mtp`` divides by a density ratio.
#: ``att`` and ``atc`` also read ``g_bounds_conditional`` where it applies, which is a
#: second bound this row does not quote. A fit that targets none of the groups below
#: therefore reports its factor rows and no product.
_COMPOSED_GROUPS = frozenset({"mean", "regime", "msm"})
#: What each excluded group divides by instead, so that the omission explains itself. A
#: silently absent row reads exactly like a fit with nothing to report, which is the one
#: thing it must not be confused with: the reader cannot tell "your denominator is fine"
#: from "this diagnostic does not cover your estimand" unless the report says which.
_COMPOSED_EXCLUSIONS: dict[str, str] = {
"att": "P(A=a) rather than g(W), and reweights the reference arm by the propensity odds",
"atc": "P(A=a) rather than g(W), and reweights the other arm by the propensity odds",
"ipsi": "a tilted mechanism that discards the propensity",
"mtp": "a density ratio rather than an arm probability",
}
[docs]
@dataclass(frozen=True)
class PositivityReport:
"""Overlap diagnostics for a fitted TMLE.
Parameters
----------
propensity_quantiles : dict of str to dict of float to float
Quantiles of ``g(W)``, overall and within each treatment arm.
tail_mass : dict of float to dict of str to float
Fraction of units with ``g(W)`` below each threshold, and above its mirror.
effective_sample_size : dict of str to dict of str to float
Kish ESS of the inverse-probability weights per arm, with the nominal arm
size alongside for comparison.
weight_share : dict of str to dict of str to float
Share of the total weight held by the largest 1% and 5% of weights.
truncated : dict of str to float
Count and fraction of propensity scores clipped by the truncation bounds, and
the most extreme untruncated value.
clever_covariate_max : dict of str to float
Largest reconstructed absolute clever-covariate value per targeted estimand
family. This is a covariate-scale diagnostic, not a score contribution or an
observation's influence on the estimate.
group_score_load : dict of str to GroupScoreLoadRow
Absolute residual-multiplier concentration for each targeted family. A row holds
fourteen keys: ``equation``, ``n_total``, ``n_targeted``, ``effective``,
``targeted_ratio``, ``total_ratio``, ``top_1pct``, ``top_5pct``, ``max_load``,
``zero_load``, ``lower_bound``, ``upper_bound``, ``clipped_count`` and
``clipped_fraction``. It carries no draw count, which
:attr:`PositivityReport.n_repeats` records once for the whole fit. The fitted
fluctuation retains ``abs(w_i * h_ij)`` on the
exact rows and with the exact weights used to evaluate each score equation. This
report computes one Kish
concentration ratio per column, then retains the most concentrated column and
names it under ``equation``. Every metric in a group row comes from that one
column. It never sums columns with different units.
``effective`` is the Kish-equivalent row count of that absolute load. It is a
descriptive concentration index. It is not the estimator's information, actual
influence, or the sample size used for its standard error. ``targeted_ratio``
divides it by the rows in the fitted score mask. ``total_ratio`` divides it by
the full data row count, so outcome missingness and intermediate selection remain
visible. The two top shares, ``max_load``, and ``zero_load`` describe the same
selected column.
``n_targeted`` counts rows in the fitted score mask, including a row whose selected
equation has a structural zero. It does not count nonzero residual contributions.
``lower_bound``, ``upper_bound``, ``clipped_count``, and ``clipped_fraction``
describe the treatment-mechanism bound used by that group. ATT and ATC therefore
report ``g_bounds_conditional`` rather than the marginal bound.
The groups that reach this table through
:meth:`~cleverly.assessment.DiagnosticsFacade.support` are ``mean``, ``att``,
``atc`` and ``msm``. A regime, shift or incremental fit gets its own support
report from that method instead and so gets no row here; ``regime`` and ``ipsi``
appear only when :func:`positivity_report` is called directly. ``mtp`` never
appears at all, because ``policies=`` needs a continuous treatment and this report
refuses one.
group_score_load_omissions : dict of str to str
Machine-readable reasons an exact group diagnostic or one of its components could
not be reported. Exact score-weight loads are never reconstructed from a generic
submodel, and clipping for a targeted treatment mechanism is unavailable unless
its fitted clipping mask was retained.
bounds : tuple of float
Truncation bounds the fit applied to the treatment mechanism.
n : int
Number of observations the diagnostics were computed over.
mechanisms : dict of str to dict of str to float
Overlap for the *other* denominators in the clever covariate:
``P(Delta = 1 | A, W)`` when outcomes are missing, and ``P(Z = z | A, W)`` for a
controlled direct effect. Every row carries the same nine keys: ``min``, ``q01``,
``q05`` and ``median`` of the mechanism as fitted; ``clipped`` and
``clipped_fraction`` for the cells the bound moved; and ``ess_ratio``,
``top_1pct`` and ``top_5pct`` for the load the weights leave behind. Empty when
neither applies.
The three load measures are taken over the weights the estimating equation
actually forms: the *bounded* mechanism at each unit's realised arm, times that
unit's observation weight, over the rows whose residual it multiplies. They are
on the same scale as the propensity's :attr:`effective_sample_size` and
:attr:`weight_share`, and fold in the design for the same reason those do, so the
two tables can be read side by side.
These deserve reporting for exactly the reason ``g`` does: they enter the
estimating equation as a denominator, so a value near zero can give one
observation a large residual multiplier. Unlike ``g`` they are one-sided -- only
the approach to zero matters -- and they are easy to overlook, because a fit can
have perfectly
healthy propensity overlap and still be resting on a handful of rows that were
very unlikely to be observed at all.
**A derived row reports the complete denominator** for a missing-outcome fit
(``P(A=a,Delta=1|W)``), a controlled direct effect (``P(A=a,Z=z|W)``), or both at
once (``P(A=a,Delta=1,Z=z|W)``). The ``z`` in those names is the level the fit
targets, which the factor row beside them states outright. It is not a further
fitted mechanism and it was never held to a bound of its own: the estimator
truncates each factor separately and multiplies. So its ``min`` and quantiles
stay on the untruncated product, as the fitted rows' do, while its ``clipped``
counts the cells where *any* factor moved -- which is the truncation this
denominator underwent, and which a comparison of the two products would miss
wherever a zero in one factor hides a clip in another. A per-arm outcome-adaptive
fit regresses each factor on ``Qbar(a, W)``, so its row names that design as the
condition: ``P(A=a,Delta=1|Qbar(a,W))``, with ``,S`` added under baseline strata.
The derived row is the reason this table exists. Each factor can look
comfortable on its own while their product concentrates the weight, because the
denominator factors multiply. That does not imply their Kish ratios multiply.
Only the product is what the clever covariate divides by.
It describes the fit's **marginal-mean** estimands -- the ``ey`` family, its
contrasts, a regime, and a marginal structural model. A fit that targets none of
those gets its factor rows and no product row, because an ATT or ATC covariate
divides by ``P(A = a)`` rather than ``g_a(W)`` and an incremental one drops the
propensity entirely.
composed_excluded : tuple of str
Targeted groups the derived row does not describe, on a fit that has the factors
to build one. Empty when every targeted group forms the product, and empty when
no factor stands beside ``g`` at all, because then there is no derived row for any
group to be outside of.
This is here so the omission explains itself. A silently absent row reads exactly
like a fit with nothing to report, and those are the two readings that must never
be confused: one says the denominator is healthy, the other says this diagnostic
does not cover the estimand asked for. :meth:`summary` names these groups and
what each divides by instead.
nuisance_bound : float
The lower bound applied to the *fitted* mechanisms above. A derived row has no
single such bound -- see :meth:`_bound_label`.
simplex_deviation : float
Largest ``|sum_a g(a | W) - 1|`` across rows *after* truncation, and ``0`` for a
two-armed fit, where the complement form preserves the sum exactly.
Non-zero is expected rather than alarming: with more than two arms the bounds are
applied arm by arm and deliberately **not** renormalised, because rescaling a row
back onto the simplex can push a column below the floor and so undo the only
thing truncation is for. The number is reported because it is the size of that
deliberate inconsistency, and a large value says the bounds are binding hard --
which is a positivity finding, not a bookkeeping one. It does not bias the
plug-in: the plug-in averages targeted predictions and contains no mechanism at
all.
n_repeats : int
Cross-fitting draws the fit combined. Everything else here describes
the **first** of them, because overlap is a property of one fitted
mechanism and a combination of draws would describe one no estimate came from.
backend : str or None
Dataframe backend :meth:`to_frame` returns when ``data`` is omitted.
treatment_mechanism : {"estimated", "known"}
Whether ``g(W)`` is an estimate or the mechanism the data declares known. A known
mechanism never moves under the fitted bounds, because the fit refuses a bound pair
that would move it.
"""
propensity_quantiles: dict[str, dict[float, float]]
tail_mass: dict[float, dict[str, float]]
effective_sample_size: dict[str, dict[str, float]]
weight_share: dict[str, dict[str, float]]
truncated: dict[str, float]
clever_covariate_max: dict[str, float]
bounds: tuple[float, float]
n: int
mechanisms: dict[str, dict[str, float]] = field(default_factory=dict)
composed_excluded: tuple[str, ...] = ()
nuisance_bound: float = 0.0
simplex_deviation: float = 0.0
#: How many cross-fitting draws the fit combined. Everything above describes the
#: **first** of them, and this is here so a reader knows that. Overlap is a property
#: of one fitted mechanism, and combining ``R`` propensity vectors would produce a
#: perfectly good estimate of ``g`` that is nonetheless not the object any reported
#: ``psi`` was computed from -- a different aggregation from the one the estimates use,
#: under the same heading. The draws share the data and differ only in the split, so
#: their overlap is near identical in practice; when it is not, that is itself worth
#: seeing rather than averaging away.
n_repeats: int = 1
#: 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
group_score_load: dict[str, GroupScoreLoadRow] = field(default_factory=dict)
#: A fluctuation without exact per-column absolute score weights is omitted rather
#: than reconstructed from a submodel that might not be the one the fit used (for
#: example under fold-specific targeting).
group_score_load_omissions: dict[str, str] = field(default_factory=dict)
treatment_mechanism: str = "estimated"
[docs]
def to_frame(self, data: Any = None) -> Any:
"""Propensity quantiles as a tidy frame.
Parameters
----------
data : Any
A dataframe or fitted container whose backend to match. ``None`` uses
:attr:`backend`.
Returns
-------
dataframe
One row per ``(arm, quantile)`` of the treatment mechanism.
"""
rows: list[tuple[str, float, float]] = [
(group, quantile, value)
for group, quantiles in self.propensity_quantiles.items()
for quantile, value in quantiles.items()
]
payload = {
"group": [row[0] for row in rows],
"quantile": [row[1] for row in rows],
"propensity": [row[2] for row in rows],
}
return emit_frame(payload, data, backend=self.backend)
def _bound_label(self, name: str) -> str:
"""How a mechanism row was truncated, as text a message can print.
A derived row is the only kind with **more than one** bound: the estimator clips
each factor separately and multiplies them, so naming any one alone -- which both
call sites used to do, always ``nuisance_bound`` -- quotes a bound that row was
never held to. How many it names comes from
:data:`_COMPOSED_NUISANCE_FACTORS`, so a further composition does not need a
second place to agree with.
"""
if name in _COMPOSED_MECHANISMS:
nuisance_factors = _COMPOSED_NUISANCE_FACTORS[name]
tail = " x ".join(f"[{self.nuisance_bound:.4g}, 1]" for _ in range(nuisance_factors))
return f"[{self.bounds[0]:.4g}, {self.bounds[1]:.4g}] x {tail}, factor by factor"
return f"[{self.nuisance_bound:.4g}, 1]"
[docs]
def summary(self) -> str:
"""A printable overlap report.
Returns
-------
str
A printable report, one line per reported quantity.
"""
lines = [
"Positivity / overlap diagnostics",
"-" * 32,
f"n = {self.n}; propensity truncated to [{self.bounds[0]:.4g}, {self.bounds[1]:.4g}]"
+ (
"; treatment mechanism known (declared)"
if self.treatment_mechanism == "known"
else ""
),
]
if self.n_repeats > 1:
lines.append(
f"describing {format_draw(REPORTED_DRAW, self.n_repeats)}: overlap is a "
"property of one fitted mechanism, not of the median-combined estimate"
)
lines.append("")
quantiles = sorted(next(iter(self.propensity_quantiles.values())))
lines.append(
format_table(
["g(W) quantile", *[f"{q:.0%}" for q in quantiles]],
[
[group, *[f"{values[q]:.4f}" for q in quantiles]]
for group, values in self.propensity_quantiles.items()
],
)
)
lines.append("")
lines.append(
format_table(
["arm", "n", "effective n", "ESS / n", "top 1% weight", "top 5% weight"],
[
[
arm,
f"{ess['n']:.0f}",
f"{ess['effective']:.1f}",
f"{ess['ratio']:.3f}",
f"{self.weight_share[arm]['top_1pct']:.3f}",
f"{self.weight_share[arm]['top_5pct']:.3f}",
]
for arm, ess in self.effective_sample_size.items()
],
)
)
lines.append("")
lines.append(
format_table(
["threshold", "P(g < t)", "P(g > 1-t)"],
[
[f"{threshold:.3g}", f"{mass['below']:.4f}", f"{mass['above']:.4f}"]
for threshold, mass in sorted(self.tail_mass.items())
],
)
)
lines.append("")
lines.append(
f"truncated: {self.truncated['count']:.0f} unit(s) "
f"({self.truncated['fraction']:.2%}); most extreme untruncated g(W) = "
f"{self.truncated['most_extreme']:.5g}"
)
if self.group_score_load:
lines.append("")
lines.append(
format_table(
[
"group",
"equation",
"target rows",
"Kish-equivalent rows",
"Kish / target",
"Kish / total",
"top 1% load",
"top 5% load",
"max |w h|",
"zero load",
"g bound",
"clipped",
],
[
[
group,
str(load["equation"]),
f"{load['n_targeted']:.0f}/{load['n_total']:.0f}",
f"{load['effective']:.1f}",
f"{load['targeted_ratio']:.3f}",
f"{load['total_ratio']:.3f}",
f"{load['top_1pct']:.3f}",
f"{load['top_5pct']:.3f}",
f"{load['max_load']:.4g}",
f"{load['zero_load']:.0f}",
(
f"[{load['lower_bound']:.4g}, {load['upper_bound']:.4g}]"
if np.isfinite(load["lower_bound"])
else "not applicable"
),
(
f"{load['clipped_count']:.0f} ({load['clipped_fraction']:.2%})"
if np.isfinite(load["clipped_fraction"])
else "unavailable"
),
]
for group, load in self.group_score_load.items()
],
)
)
lines.append(
"(absolute residual-multiplier load from the fitted equation; Kish values describe "
"residual-multiplier concentration, not residual contributions or estimator "
"information, influence, or precision; target rows are score-mask rows and "
"include structural zeros)"
)
if self.clever_covariate_max:
for group, value in self.clever_covariate_max.items():
lines.append(f"max |clever covariate| ({group}): {value:.4g}")
if self.group_score_load_omissions:
lines.append("")
for group, reason in self.group_score_load_omissions.items():
lines.append(f"group diagnostic omitted ({group}): {reason}")
if self.mechanisms:
lines.append("")
lines.append(
format_table(
[
"mechanism",
"min",
"1%",
"5%",
"median",
"ESS / n",
"top 1% weight",
"top 5% weight",
"clipped",
],
[
[
name,
f"{stats['min']:.4f}",
f"{stats['q01']:.4f}",
f"{stats['q05']:.4f}",
f"{stats['median']:.4f}",
f"{stats['ess_ratio']:.3f}",
f"{stats['top_1pct']:.3f}",
f"{stats['top_5pct']:.3f}",
f"{stats['clipped']:.0f} ({stats['clipped_fraction']:.2%})",
]
for name, stats in self.mechanisms.items()
],
)
)
truncations = "; ".join(
f"{name} truncated to {self._bound_label(name)}" for name in self.mechanisms
)
lines.append(f"({truncations}; each row counts both arms)")
for name in sorted(_COMPOSITE_ROWS & self.mechanisms.keys()):
stats = self.mechanisms[name]
lines.append(
f"{name}: minimum {stats['min']:.4g}; "
f"{stats['below_g_lower_fraction']:.2%} of unit-arm cells below the "
f"propensity floor {self.bounds[0]:.4g}. The composite is tilted inside "
f"[{stats['composite_lower']:.4g}, {self.bounds[1]:.4g}], the product of "
"the factor floors."
)
note = self._composed_note()
if note is not None:
lines.append("")
lines.append(note)
lines.append("")
lines.append(self.verdict())
return "\n".join(lines)
def _composed_note(self) -> str | None:
"""The sentence a fit gets when a derived row does not cover what it targeted.
Returns ``None`` when there is nothing to say: every targeted group forms the
product, or no fitted factor stands beside ``g`` so there is no derived row at
all. Otherwise it names the groups and what each divides by instead, because an
omission that does not explain itself reads exactly like a clean bill of health.
The refusal stands, and the note now says where those groups *are* covered. The
per-group score-load table reads each covariate as built rather than a denominator
it never forms, so an excluded group has a load measure even though it has no
derived row. Naming it keeps the sentence from reading as "nothing here for you".
"""
if not self.composed_excluded:
return None
excluded = "; ".join(
f"{group}, which divides by "
+ _COMPOSED_EXCLUSIONS.get(group, "a denominator this row does not form")
for group in self.composed_excluded
)
composed = sorted(_COMPOSED_MECHANISMS & self.mechanisms.keys())
lead = (
f"{composed[0]} does not describe"
if composed
else "no derived denominator row is reported for"
)
pointer = (
" Their score-weight load is in the group table above." if self.group_score_load else ""
)
return f"({lead} {excluded}. Read the factor rows above for those estimands.{pointer})"
@property
def severity(self) -> Literal["adequate", "strain", "serious"]:
"""Which tier this report's own thresholds place the fit in.
Returns
-------
{"adequate", "strain", "serious"}
The tier behind :meth:`verdict`, for a caller that must branch on it.
"""
return self._verdict_parts()[0]
[docs]
def verdict(self) -> str:
"""A one-line reading of the diagnostics.
Returns
-------
str
One line saying whether the overlap supports the reported estimate.
"""
return self._verdict_parts()[1]
def _verdict_parts(self) -> tuple[Literal["adequate", "strain", "serious"], str]:
"""Decide the tier and the sentence that reports it, in one place.
A reader reads :meth:`verdict`; the combined assessment row reads
:attr:`severity`. Both come from here because the alternative is what shipped:
the aggregate row re-derived its own thresholds and so reported ``passed`` for a
fit this method called a serious positivity problem. No caller outside this
method may restate the numbers below.
**The effective-sample-size ratio is reported and never graded.** A Kish ratio
is a descriptive quantity with no derived cutoff, so a tier keyed on it would be
a house convention presented as a finding, and a reader would take ``passed`` as
a positivity clearance this package cannot give. Every branch states the share
instead, and the analyst judges it against the question. What is graded is the
truncated fraction, which is not a judgement about the data: it counts the rows
the fit clipped at a bound the caller configured, and therefore identifies where
the score used a bounded rather than fitted mechanism value.
The **group** quantity from :attr:`group_score_load` is a different descriptive
object: concentration of ``abs(w_i * h_ij)`` in one fitted score equation. It is
stated separately, never pooled with arm or mechanism ESS, and adds no tier or
threshold. In particular it is not estimator effective sample size, information,
actual influence, precision, or itself a positivity test.
"""
# Non-finite ratios are dropped rather than compared. An arm with no rows stores
# NaN here by construction, and `min` over a sequence containing NaN returns
# whichever value it met first, so an unfiltered comparison made this report's
# reading depend on arm order.
ratios = [
float(ess["ratio"])
for ess in self.effective_sample_size.values()
if np.isfinite(ess["ratio"])
]
fraction = float(self.truncated["fraction"])
clip_label = (
f"the marginal treatment-mechanism bound [{self.bounds[0]:.4g}, {self.bounds[1]:.4g}]"
)
for group, load in self.group_score_load.items():
group_fraction = float(load["clipped_fraction"])
if np.isfinite(group_fraction) and group_fraction > fraction:
fraction = group_fraction
clip_label = (
f"group {group!r}'s treatment-mechanism bound "
f"[{load['lower_bound']:.4g}, {load['upper_bound']:.4g}]"
)
share = (
(
f"The narrowest arm has a Kish-equivalent weight count of {min(ratios):.0%} "
"of its nominal rows. This is concentration, not estimator information, and "
"no threshold is applied because none is derived."
)
if ratios
else "No arm reports a finite effective sample size."
)
group_rows = [
(group, load)
for group, load in self.group_score_load.items()
if np.isfinite(load["targeted_ratio"])
]
if group_rows:
group, load = min(group_rows, key=lambda item: item[1]["targeted_ratio"])
share += (
f" Absolute-load concentration is greatest for group {group!r}, equation "
f"{load['equation']!r}: {load['effective']:.1f} Kish-equivalent rows out of "
f"{load['n_targeted']:.0f} score-mask rows "
f"({load['targeted_ratio']:{SCORE_LOAD_VERDICT_RATIO_FORMAT}}; "
f"{load['total_ratio']:{SCORE_LOAD_VERDICT_RATIO_FORMAT}} of "
f"{load['n_total']:.0f} total rows). This is a "
"descriptive residual-multiplier concentration index. It excludes the residual "
"and is not estimator effective sample size, information, actual influence, "
"precision, or a positivity test."
)
for name, stats in self.mechanisms.items():
# Checked before the propensity verdict, because this is the failure a reader
# is least likely to be looking for: overlap in `g` can be immaculate while
# the estimate rests on a few rows that were very unlikely to be observed.
# Two triggers, one of which only reports. A clipped fraction is truncation
# and grades like the propensity's. A low mechanism effective sample size
# earns the *sentence*, because a reader who never sees the number cannot
# weigh it, but it does not move the tier -- grading it would be the invented
# cutoff this report refuses to apply to the propensity.
composed = name in _COMPOSED_MECHANISMS
# A derived row's clipped cells are the union over its factors, and its
# propensity factor is the one the branch below already reports -- accurately,
# and against the curve that can move it. Triggering here on that union put a
# mechanism's name on a verdict whose every number came from `g`, and sent the
# reader to `mechanism=True`, which sweeps `nuisance_bound` alone. So a derived
# row earns its sentence on joint leverage only, which is the one thing neither
# its factor rows nor the propensity branch can show. Each factor's own
# clipping still triggers on that factor's own row.
clipped = not composed and stats["clipped_fraction"] > 0.01
if clipped or stats["ess_ratio"] < 0.6:
leverage = (
"It is the whole denominator of the clever covariate for this fit's "
"marginal-mean estimands, so those rows carry outsized leverage "
"however each factor looks on its own."
if composed
else "It divides the clever covariate exactly as g(W) does, so those "
"rows carry outsized leverage whatever the propensity overlap looks like."
)
# Either bound can have moved a derived row's cells, so naming one curve
# would send the reader to a knob that cannot move what they were shown.
advice = (
"Check truncation_curve() and truncation_curve(mechanism=True); this "
"row is clipped by either bound."
if composed
else "Check truncation_curve(mechanism=True)."
)
return (
"serious"
if fraction > 0.05
else "strain"
if fraction > 0.01 or clipped
else "adequate",
f"VERDICT: {name} strains the estimate. It falls to {stats['min']:.4g} at "
f"its smallest and leaves an effective {stats['ess_ratio']:.0%} of the "
f"rows it weights ({stats['clipped_fraction']:.2%} clipped at "
# `share` is carried onto this branch too, so that "every verdict
# states the arm share and the group share" holds on the tier that
# is already flagging strain rather than only on the three below.
f"{self._bound_label(name)}). Maximum treatment-mechanism clipping is "
f"{fraction:.2%} at {clip_label}. {leverage} {advice} {share}",
)
if fraction > 0.05:
return (
"serious",
f"VERDICT: truncation is carrying this estimate. {fraction:.1%} of units were "
f"clipped at {clip_label}, so the fitted score used bounded rather than fitted "
"mechanism values for those units. Treat "
"the estimate as sensitive to this finite-sample regularisation. The bound "
"does not change the requested estimand. Check "
f"truncation_curve() before drawing conclusions. {share}",
)
if fraction > 0.01:
return (
"strain",
f"VERDICT: some truncation. {fraction:.1%} of units were clipped at "
f"{clip_label}. "
"Report truncation_curve() alongside the estimate so readers can see how much "
f"the answer depends on it. {share}",
)
return (
"adequate",
f"VERDICT: no truncation-driven fragility detected ({fraction:.1%} of units "
f"clipped at {clip_label}). {share}",
)
def __repr__(self) -> str: # pragma: no cover - debugging aid
return self.summary()
[docs]
def positivity_report(result: TMLEResult) -> PositivityReport:
"""Compute overlap diagnostics for a fitted TMLE.
Two arms and more than two are reported by separate functions rather than one
parameterised by arm count, because the *questions* differ. With two arms there is
one propensity and overlap is symmetric: ``g`` near zero and ``g`` near one are the
same problem seen from either arm, and the interesting split is treated versus
control. With more there is no single margin and no mirror -- each arm has its own
denominator, which has to be reported and truncated in its own right. Collapsing
the two into one function would mean picking definitions that read oddly in both.
A **continuous** treatment is refused rather than given a third branch. Every field
of :class:`PositivityReport` is per arm -- quantiles of ``g(a | W)``, tail mass,
effective sample size, weight share -- and a dose has none, so the report would come
back empty with a ``simplex_deviation`` of ``1.0`` computed from a zero-column
mechanism: the largest value the field can take, reported as a finding. The question
a shift fit actually has to answer is about the *density ratio* at the shifted dose,
which :func:`~cleverly.interventions.check_policy_support` answers.
Parameters
----------
result : TMLEResult
A fitted result.
Returns
-------
PositivityReport
Overlap read off the stored nuisance fits, without refitting.
Raises
------
CapabilityError
If the fit estimated no treatment mechanism. A missing-outcome
``NaturalCourseMean`` fit has no propensity for an arm-wise support report to
read; inspect the missingness row of ``diagnostics.nuisance_models()`` instead.
DataError
If the treatment is continuous, which has no per-arm propensity either.
"""
if not result.nuisance.fits_treatment:
raise CapabilityError.from_refusal(NATURAL_COURSE_SUPPORT_REFUSAL)
if result.data.is_continuous_treatment:
raise DataError(
f"{result.data.treatment_name} is continuous, so there is no per-arm "
"propensity to tabulate and this report has no rows to fill. A shift's "
"positivity question is whether the density ratio g(a - delta | W) / "
"g(a | W) stays bounded, not whether an arm probability does.",
reason=RefusalReason.DIFFERENT_ESTIMAND,
alternative=(
"check_policy_support answers it, which res.diagnostics.support() reaches "
"for a fit that declared policies=. Reaching this report instead means none "
"were declared. With delta= or intermediate= there is a mechanism in the "
"denominator as well, and it is the one bound this axis actually has: "
"res.diagnostics.truncation_curve(mechanism=True) sweeps it."
),
)
if result.data.is_binary_treatment:
return _binary_positivity_report(result)
return _multi_arm_positivity_report(result)
def _binary_positivity_report(result: TMLEResult) -> PositivityReport:
"""Overlap for a two-armed treatment, in terms of the single propensity ``g(W)``."""
data = result.data
bounds = result.config.g_bounds
propensity = result.nuisance.propensity
raw = propensity.arm(1.0)
# Each arm is the rows that recorded it. A row with a missing treatment is in neither,
# so the control arm is ``A == 0`` rather than the complement of the treated rows.
treated = data.treatment == 1.0
control = data.treatment == 0.0
quantiles: dict[str, dict[float, float]] = {
"overall": {q: float(np.quantile(raw, q)) for q in _QUANTILES},
"treated": {q: float(np.quantile(raw[treated], q)) for q in _QUANTILES},
"control": {q: float(np.quantile(raw[control], q)) for q in _QUANTILES},
}
tail_mass = {
threshold: {
"below": float(np.mean(raw < threshold)),
"above": float(np.mean(raw > 1.0 - threshold)),
}
for threshold in _THRESHOLDS
}
# Both denominators come out of the one truncation the estimator applies, rather than
# from a clip written here and a complement written beside it. The two forms agree
# only while the bound pair is symmetric.
truncation = propensity.truncate(bounds)
bounded = truncation.values
ess: dict[str, dict[str, float]] = {}
share: dict[str, dict[str, float]] = {}
for arm, mask, weights in (
("treated", treated, 1.0 / bounded[:, propensity.column_for(1.0)]),
("control", control, 1.0 / bounded[:, propensity.column_for(0.0)]),
):
arm_weights = weights[mask] * data.weights[mask]
ess[arm] = {
"n": float(mask.sum()),
"effective": _kish_ess(arm_weights),
"ratio": _kish_ess(arm_weights) / float(mask.sum()) if mask.any() else float("nan"),
}
share[arm] = {
"top_1pct": top_weight_share(arm_weights, 0.01),
"top_5pct": top_weight_share(arm_weights, 0.05),
}
clipped = truncation.units
if propensity.simplex:
# On the simplex the control denominator is the complement of ``g1``, so the most
# extreme value is the nearer of ``g1`` to zero and to one.
inside = raw[~clipped]
most_extreme = float(min(inside.min(), 1.0 - inside.max())) if inside.size else float("nan")
else:
# Off the simplex each column is its own denominator and is clipped on its own, so
# the most extreme value is the smallest untruncated cell of either column.
cells = np.asarray(propensity.values, dtype=float)[~truncation.clipped]
most_extreme = float(cells.min()) if cells.size else float("nan")
group_score_load, group_omissions = _group_score_load(result)
return PositivityReport(
propensity_quantiles=quantiles,
tail_mass=tail_mass,
effective_sample_size=ess,
weight_share=share,
truncated={
"count": float(clipped.sum()),
"fraction": float(clipped.mean()),
"most_extreme": most_extreme,
},
clever_covariate_max={
group: _max_abs_covariate(result, group) for group in result.fluctuations
},
bounds=bounds,
n=data.n,
mechanisms=_mechanism_overlap(result),
composed_excluded=_composed_excluded(result),
nuisance_bound=result.config.missingness_bound,
n_repeats=result.n_repeats,
backend=data.backend,
group_score_load=group_score_load,
group_score_load_omissions=group_omissions,
treatment_mechanism=result.nuisance.treatment_mechanism,
)
def _multi_arm_positivity_report(result: TMLEResult) -> PositivityReport:
r"""Overlap for a ``K``-armed treatment, arm by arm.
Each arm's probability :math:`g_a(W)` is summarised over **all** rows, not just the
rows in that arm. That is the distribution positivity actually depends on: the
clever covariate divides by :math:`g_a` and the plug-in evaluates
:math:`\bar Q(a, W)` at every unit, so a unit that could never have received arm
``a`` is a problem for arm ``a`` whichever arm it did receive.
``tail_mass`` loses its mirror here. With two arms ``g > 1 - t`` is the same
statement as ``g < t`` read from the control arm; with more arms an arm being
*likely* is not another arm being unlikely, so ``below`` counts unit-arm pairs whose
probability sits under the threshold and ``above`` counts the unit-arm pairs that are
nearly deterministic -- related, but no longer the same number counted twice.
"""
data = result.data
bounds = result.config.g_bounds
propensity = result.nuisance.propensity
raw = np.asarray(propensity.values, dtype=float)
labels = {arm: str(data.arm_label(arm)) for arm in propensity.arms}
quantiles: dict[str, dict[float, float]] = {
f"g[{labels[arm]}]": {q: float(np.quantile(propensity.arm(arm), q)) for q in _QUANTILES}
for arm in propensity.arms
}
tail_mass = {
threshold: {
"below": float(np.mean(raw < threshold)),
"above": float(np.mean(raw > 1.0 - threshold)),
}
for threshold in _THRESHOLDS
}
truncation = propensity.truncate(bounds)
bounded = truncation.values
ess: dict[str, dict[str, float]] = {}
share: dict[str, dict[str, float]] = {}
for arm in propensity.arms:
mask = data.treatment == arm
column = bounded[:, propensity.column_for(arm)]
arm_weights = (1.0 / column)[mask] * data.weights[mask]
nominal = float(mask.sum())
ess[labels[arm]] = {
"n": nominal,
"effective": _kish_ess(arm_weights),
"ratio": _kish_ess(arm_weights) / nominal if mask.any() else float("nan"),
}
share[labels[arm]] = {
"top_1pct": top_weight_share(arm_weights, 0.01),
"top_5pct": top_weight_share(arm_weights, 0.05),
}
clipped_cells = truncation.clipped
clipped_units = truncation.units
inside = raw[~clipped_cells]
# Only the approach to zero matters per arm, so the "most extreme" untruncated value
# is the smallest surviving probability rather than the two-sided minimum the binary
# report uses.
most_extreme = float(inside.min()) if inside.size else float("nan")
group_score_load, group_omissions = _group_score_load(result)
return PositivityReport(
propensity_quantiles=quantiles,
tail_mass=tail_mass,
effective_sample_size=ess,
weight_share=share,
truncated={
"count": float(clipped_units.sum()),
"fraction": float(clipped_units.mean()),
"most_extreme": most_extreme,
},
clever_covariate_max={
group: _max_abs_covariate(result, group) for group in result.fluctuations
},
bounds=bounds,
n=data.n,
mechanisms=_mechanism_overlap(result),
composed_excluded=_composed_excluded(result),
nuisance_bound=result.config.missingness_bound,
simplex_deviation=float(np.max(np.abs(bounded.sum(axis=1) - 1.0))),
n_repeats=result.n_repeats,
backend=data.backend,
group_score_load=group_score_load,
group_score_load_omissions=group_omissions,
treatment_mechanism=result.nuisance.treatment_mechanism,
)
def _composed_excluded(result: TMLEResult) -> tuple[str, ...]:
"""Targeted groups a derived row would not describe, on a fit that could build one.
The mirror of the gate in :func:`_mechanism_overlap`, and it reads the same two
conditions: a fit with no fitted factor beside ``g`` has no derived row to be outside
of, and a group in :data:`_COMPOSED_GROUPS` is described by the row rather than
excluded from it. Everything left over is what the report has to account for.
"""
nuisance = result.nuisance
has_factor = (
nuisance.missingness is not None
or nuisance.treatment_observation is not None
or (nuisance.intermediate is not None and result.intermediate_value is not None)
)
if not has_factor:
return ()
return tuple(group for group in result.fluctuations if group not in _COMPOSED_GROUPS)
def _mechanism_overlap(result: TMLEResult) -> dict[str, dict[str, float]]:
"""Overlap for the denominators other than ``g`` -- ``pi`` and the intermediate density.
Two views, because they answer different questions. The quantiles pool every arm,
since every column is used: each arm's covariate divides by the mechanism at that
arm, so the union is the set of values that appear as denominators anywhere. The
effective sample size instead takes the weights the estimating equation *actually*
forms -- ``1 / pi`` at each unit's realised arm, over the rows whose residual it
multiplies -- and is reported on the same scale as the propensity's ESS so that the
two can be read side by side.
Beside each fitted mechanism the report carries the *product* the clever covariate
divides by, named in :data:`_COMPOSED_ROWS`. It is built here rather than stored on
the nuisance state, so a cached product cannot go stale against the factors it came
from, and each factor is bounded through the accessor that owns its rule before the
multiplication -- which is the order the estimator uses.
Which rows those are depends on the estimand. A row with no recorded outcome
contributes a genuine zero to the residual term, and so does a row whose intermediate
is not the level being targeted; neither is weighted by any mechanism, so neither
belongs in an effective sample size. For a controlled direct effect that is a real
difference rather than a technicality -- roughly half the sample is typically at the
other level, and averaging it in would report an ESS for a weighting that never
happened.
"""
data = result.data
nuisance = result.nuisance
missingness_bound = result.config.missingness_bound
# A row with an unrecorded treatment contributes a genuine zero to the residual term,
# as a row with an unrecorded outcome does.
contributing = (
np.asarray(targeted_rows(data, result.intermediate_value), dtype=bool)
& data.treatment_recorded
)
out: dict[str, dict[str, float]] = {}
# Which column of an `(n, K)` mechanism each row realised. Written as a gather rather
# than as `np.where(treatment == 1, values[:, 1], values[:, 0])`, which reports arm 0's
# denominator for every unit outside arms 0 and 1 -- silently, and on a report that is
# otherwise per arm. A two-armed fit gathers the same two columns that expression did.
realised = np.zeros(data.n, dtype=int)
for arm in nuisance.propensity.arms:
realised[data.treatment == arm] = nuisance.propensity.column_for(arm)
rows = np.arange(data.n)
def summarize(
raw: FloatArray, bounded: FloatArray, clipped: BoolArray | None = None
) -> dict[str, float]:
"""One mechanism row: raw quantiles, bounded-weight load, and the cells that moved.
The two arrays are deliberately different. Quantiles pool every column, because
every column appears as a denominator somewhere, and they describe the mechanism
as fitted. The weights are the ones the estimating equation actually forms: the
*bounded* array at each unit's realised arm, times that unit's observation weight,
over the rows whose residual it multiplies. ``clipped`` defaults to the cells the
bound moved, which is the rule for a single factor; a derived row passes the union
of its factors' masks instead, because its own product was never bounded.
"""
flat = raw.reshape(-1)
# A factor that does not depend on the arm is one column, read at every row.
at_arm = bounded if bounded.ndim == 1 else bounded[rows, realised]
weights = data.weights[contributing] / at_arm[contributing]
changed = raw != bounded if clipped is None else clipped
return {
"min": float(flat.min()),
"q01": float(np.quantile(flat, 0.01)),
"q05": float(np.quantile(flat, 0.05)),
"median": float(np.median(flat)),
"ess_ratio": (
_kish_ess(weights) / float(weights.size) if weights.size else float("nan")
),
"top_1pct": top_weight_share(weights, 0.01),
"top_5pct": top_weight_share(weights, 0.05),
"clipped": float(np.count_nonzero(changed)),
"clipped_fraction": float(np.mean(changed)),
}
# Each fitted nuisance factor, as the pair the report needs: the array as fitted, and
# the same array as the estimator bounds it at targeting time. Both come from the
# accessors that own those rules -- `bounded_missingness`, and `intermediate_density`,
# which is the one place the `z` / `1 - z` convention lives, so the row reports the
# level the covariate actually divides by rather than its complement.
level = result.intermediate_value
has_observation = nuisance.missingness is not None
has_intermediate = nuisance.intermediate is not None and level is not None
has_recording = nuisance.treatment_observation is not None
outcome_label, treatment_label = _observation_labels(result)
factors: list[tuple[str, FloatArray, FloatArray]] = []
if has_recording:
recorded_raw = np.asarray(nuisance.treatment_observation, dtype=float)
factors.append(
(
treatment_label,
recorded_raw,
np.clip(recorded_raw, missingness_bound, 1.0),
)
)
if has_observation:
factors.append(
(
outcome_label,
np.asarray(nuisance.missingness, dtype=float),
np.asarray(nuisance.bounded_missingness(missingness_bound), dtype=float),
)
)
if has_intermediate and level is not None:
factors.append(
(
f"P(Z={level:.0f}|A,W)",
np.asarray(nuisance.intermediate_density(level, 0.0), dtype=float),
np.asarray(nuisance.intermediate_density(level, missingness_bound), dtype=float),
)
)
for name, raw, bounded_factor in factors:
out[name] = summarize(raw, bounded_factor)
# The product is a derived denominator, not a further fitted or targeted mechanism.
# Report it without storing it on the nuisance state so the treatment, observation and
# intermediate probabilities cannot become stale relative to a cached product. It is
# reported only for the groups that form it -- see `_COMPOSED_GROUPS`, whose comment
# names what the others divide by instead. Refusing the row there is the point: an
# ATT fit's covariate never forms `g_a * pi_a`, so a row asserting that product is its
# whole denominator would describe an estimand nobody asked for.
if factors and not _COMPOSED_GROUPS.isdisjoint(result.fluctuations):
# The estimator truncates each factor **separately** and multiplies -- see
# `build_submodel`'s `1 / (g * pi * pz)` and the missing-outcome reductions -- so
# `g_bounds[0] * missingness_bound` is a floor no code applies. Counting the cells
# beneath it reports a strict subset of the cells truncation altered, because a
# small factor beside a large one leaves the product above the product of the
# floors: at the shipped defaults that floor is around 5e-4, so the row read zero
# on fits where a third of the observation mechanism was pinned. What is reported
# instead is the union of the factors' own clipped cells, which is the truncation
# this row's denominator actually underwent -- and which still marks a cell that a
# zero in another factor would have hidden from a product-level comparison. The
# propensity contributes `Truncation`'s mask rather than a rebuilt predicate,
# because two arms on the simplex are clipped through `g1` and its complement, so
# a cell can move without ever standing below the bound.
truncation = nuisance.propensity.truncate(result.config.g_bounds)
raws = [np.asarray(nuisance.propensity.values, dtype=float)]
boundeds = [np.asarray(truncation.values, dtype=float)]
masks = [np.asarray(truncation.clipped, dtype=bool)]
width = raws[0].shape[1]
for _, raw, bounded_factor in factors:
# A factor that does not depend on the arm multiplies every arm's column.
if raw.ndim == 1:
raw = np.repeat(raw[:, None], width, axis=1)
bounded_factor = np.repeat(bounded_factor[:, None], width, axis=1)
raws.append(raw)
boundeds.append(bounded_factor)
masks.append(raw != bounded_factor)
name = _composed_row(result, has_observation, has_intermediate, has_recording)
product = np.prod(np.stack(boundeds), axis=0)
out[name] = summarize(
np.prod(np.stack(raws), axis=0),
product,
np.logical_or.reduce(np.stack(masks)),
)
if name in _COMPOSITE_ROWS:
lower = float(result.config.g_bounds[0])
factor_count = int(has_observation) + int(has_recording)
out[name]["below_g_lower_fraction"] = float(np.mean(product < lower))
out[name]["composite_lower"] = lower * float(missingness_bound) ** factor_count
return out
def _kish_ess(weights: FloatArray) -> float:
"""Kish's ESS, answering ``nan`` for an arm with no rows rather than zero.
The distinction is this module's own and is why it still has a wrapper: an arm nothing
was selected for has no effective sample size to report, while an arm whose weights
sum to zero has one and it is zero. The formula itself is
:func:`~cleverly.data.weighting.effective_sample_size`.
"""
w = np.asarray(weights, dtype=float)
if w.size == 0:
return float("nan")
return effective_sample_size(w, on_degenerate=0.0)
def _group_submodel(result: TMLEResult, group: str) -> Submodel:
"""The clever covariate one targeted family was fluctuated along.
Rebuilt from the data, the nuisance estimates and the config rather than from the
estimator, so this stays a real number on a result whose estimator is gone.
Every reading of a group's covariate comes through here. The largest absolute value
and the load measures describe the same array by construction, which they would not
if each caller rebuilt the submodel with bounds and a reference arm of its own.
"""
bounds = g_bounds_for(group, result.config.g_bounds, result.config.g_bounds_conditional)
data, nuisance = result.data, result.nuisance
if result.extra.get("missing_data") == "composite":
# A composite fit was fluctuated along C_a / g_c, built at the derivation site at the
# fit's own bounds. Rebuilt here from the raw data, the covariate would omit the
# treatment observation factor and read the codes of rows that are not C_a = 1.
state = composite_state(
data, nuisance, g_bounds=bounds, nuisance_bound=result.config.missingness_bound
)
data, nuisance, bounds = state.data, state.nuisance, state.bounds
return build_submodel(
data,
nuisance,
group,
bounds=bounds,
nuisance_bound=result.config.missingness_bound,
intermediate_value=result.intermediate_value,
# The fit's own reference, so a conditional-effect covariate is rebuilt with the
# contrasts it was targeted with rather than with the lowest arm's.
reference=result.config.reference_arm,
# A working model with a non-identity link has a covariate that reads its own
# coefficients, so "the covariate" is only defined once they are named. The ones
# the fit *reports* at are the right choice and the only defensible one: they are
# the equation this fit solved, and under fold-wise targeting they are the single
# beta a per-fold covariate has no other summary of.
msm_beta=_reported_beta(result, group),
)
def _max_abs_covariate(result: TMLEResult, group: str) -> float:
"""Largest absolute clever-covariate value for one targeted family."""
return _group_submodel(result, group).max_abs
def _group_score_load(
result: TMLEResult,
) -> tuple[dict[str, GroupScoreLoadRow], dict[str, str]]:
"""Summarise the most concentrated fitted score equation in each target group.
The fluctuation artifact is authoritative. Rebuilding a generic submodel is not:
fold-targeted ATT can use fold-specific arm fractions, and a nonlinear MSM can use a
different beta in every fold. A fluctuation without the exact absolute score
weights is therefore omitted with a reason rather than approximated.
Each score column is considered separately. Kish and top-share concentration are
invariant to a scalar change of that column's units, while an L1 sum across columns is
not. The retained group row is the column with the smallest targeted-row Kish ratio,
and every metric in that row is computed from that same column.
The guards and the ten concentration values are
:func:`~cleverly.data.weighting.validate_score_loads` and
:func:`~cleverly.data.weighting.score_load_row`, shared with the intervention reports.
Adopting them adds one refusal this path did not have: a block with more rows than the
fit had is now omitted with a reason, where before it was described and published a
``total_ratio`` above one.
"""
rows: dict[str, GroupScoreLoadRow] = {}
omissions: dict[str, str] = {}
for group, fluctuation in result.fluctuations.items():
loads, reason = validate_score_loads(
getattr(fluctuation, "absolute_score_weights", None),
len(fluctuation.names),
result.data.n,
)
if loads is None:
assert reason is not None # returned together; narrows the type
omissions[group] = reason
continue
columns = [
score_load_row(loads[:, index], name, result.data.n)
for index, name in enumerate(fluctuation.names)
]
# Stable first-column selection in an exact tie. A zero-load equation has ratio
# zero and is correctly the one the report must expose rather than average away.
selected = int(np.argmin([column["targeted_ratio"] for column in columns]))
if group == "ipsi":
lower = upper = float("nan")
clipped_count = clipped_fraction = 0.0
else:
lower, upper = g_bounds_for(
group, result.config.g_bounds, result.config.g_bounds_conditional
)
if fluctuation.reduction is not None or fluctuation.mechanism is not None:
# These fits solve the outcome equation against a treatment mechanism
# changed during targeting. Its score loads are exact, but the retained
# initial nuisance cannot recover which cells the bound moved. Reporting
# initial-g clipping here would attach the right numerical bound to the
# wrong mechanism.
clipped_count = clipped_fraction = float("nan")
omissions[f"{group}.clipping"] = (
"the score used a targeted treatment mechanism whose exact clipping "
"mask was not retained"
)
else:
clipped = result.nuisance.propensity.truncate((lower, upper)).units
clipped_count = float(np.count_nonzero(clipped))
clipped_fraction = float(np.mean(clipped))
rows[group] = {
**columns[selected],
"lower_bound": float(lower),
"upper_bound": float(upper),
"clipped_count": clipped_count,
"clipped_fraction": clipped_fraction,
}
return rows, omissions
def _reported_beta(result: TMLEResult, group: str) -> Any:
"""The working model's coefficients, or ``None`` where the covariate does not read them."""
projection = getattr(result.fluctuations.get(group), "projection", None)
return None if projection is None else projection.beta
def truncation_axis(result: TMLEResult, mechanism: bool | None) -> bool:
"""Resolve the axis a truncation curve sweeps.
``None`` resolves to the only axis the fit has when it has one: the observation
mechanism on a fit that estimates no treatment law, which is the missing-outcome
natural-course mean, and ``g(W)`` on every other fit. It asks the nuisances rather
than the estimand name, because a module-level caller can hold a result that an
estimand-name test would not recognise. An explicit value is kept, so a caller who
names an axis the fit does not have hears a refusal rather than a redirection.
Parameters
----------
result : TMLEResult
A fitted point-treatment result.
mechanism : bool or None
The ``mechanism=`` argument as the caller passed it.
Returns
-------
bool
``True`` to sweep the observation or intermediate mechanism, ``False`` to sweep
``g(W)``.
"""
if mechanism is None:
return not result.nuisance.fits_treatment
return bool(mechanism)
def _refuse_observation_axis(result: TMLEResult, axis: bool) -> Refusal | None:
"""Refuse the mechanism axis on a fit with no mechanism in its clever covariate."""
fitted = (
result.nuisance.missingness is not None
or result.nuisance.intermediate is not None
or result.nuisance.treatment_observation is not None
)
if not axis or fitted:
return None
return Refusal(
RefusalReason.INVALID_DECLARATION,
"mechanism=True needs a fit with missing outcomes or an intermediate "
"variable; without one there is no mechanism in the clever covariate to "
f"truncate. Declare the indicator with {MISSING_OUTCOME_DECLARATION}, or the "
"intermediate variable with intermediate=<column> on PointTreatment or fit().",
)
def _refuse_incremental(result: TMLEResult, axis: bool) -> Refusal | None:
"""Refuse the propensity axis of an incremental fit, whose estimand contains ``g``."""
if axis or result.nuisance.incremental is None:
return None
return Refusal(
RefusalReason.DIFFERENT_ESTIMAND,
"the propensity g is *inside* the estimand for an incremental intervention, "
"so a propensity-bound curve would compare different parameters.",
alternative=(
"Use diagnostics.support(), or pass mechanism=True when a separate observation "
"mechanism was fitted."
),
)
def _refuse_natural_course(result: TMLEResult, axis: bool) -> Refusal | None:
"""Refuse the propensity axis of a fit that estimates no treatment mechanism.
Refused here rather than left to ``Propensity.truncate``, which would raise a bare
``ValueError`` naming an internal staging invariant, and only after the whole grid had
been retargeted once. The sentence names the axis the fit does have. It does not
reuse the support report's sentence, which sends the reader to ``nuisance_models()``
instead, because this caller can have the curve by changing one argument.
"""
if axis or result.nuisance.fits_treatment:
return None
return Refusal(
RefusalReason.API_REDIRECTION,
"NaturalCourseMean with missing outcomes fits no treatment propensity, so "
"there is no g(W) bound to sweep.",
alternative=(
"Pass mechanism=True to sweep the bound on P(Delta = 1 | A, W), which is the "
"only mechanism this fit truncates."
),
)
#: Every refusal of a point-treatment truncation curve that the requested axis decides, in
#: the one order the module call and the capability row both use. Each rule reads the
#: resolved axis from :func:`truncation_axis`. The names are the introspection contract,
#: so a test pins the order without respelling a message.
_TRUNCATION_RULES: tuple[tuple[str, Callable[[TMLEResult, bool], Refusal | None]], ...] = (
("observation_axis", _refuse_observation_axis),
("incremental", _refuse_incremental),
("natural_course", _refuse_natural_course),
)
def refits_reduced_regressions(result: TMLEResult) -> bool:
"""Whether a truncation curve of ``result`` refits the reduced regressions at each bound.
A guarded DR-TMLE fit alternates its targeting against the reduced regressions, and
:func:`truncation_curve` refits them inside ``retarget``. The capability row prices
that curve as a refit, and :func:`truncation_refusal` asks the estimator whether its
reduction settings admit that refit.
Parameters
----------
result : TMLEResult
A fitted point-treatment result.
Returns
-------
bool
``True`` for a DR-TMLE fit that solved its correction equations.
"""
return getattr(result, "fitted_method", None) == "drtmle" and bool(
getattr(result, "solved_corrections", False)
)
def truncation_refusal(result: TMLEResult, mechanism: bool | None = None) -> Refusal | None:
"""Return the refusal a point-treatment truncation curve meets on one axis.
:func:`truncation_curve` raises the sentence this function returns, and the
``truncation_curve`` capability row of :class:`~cleverly.assessment.DiagnosticsFacade`
quotes it for the same request. The row defers on ``mechanism`` when the default
axis is refused and the other one runs, so a combined report names the argument
instead of calling a curve the call refuses.
A guarded DR-TMLE curve retargets the cached primary nuisances and refits the reduced
regressions under the estimator's live reduction settings. It is refused on either
axis when a fit refuses those settings. A setting that only the primary split reads,
such as ``stratify_folds``, does not reach the curve and is not asked.
Parameters
----------
result : TMLEResult
A fitted point-treatment result.
mechanism : bool or None
The ``mechanism=`` argument as the caller passed it. ``None`` is resolved by
:func:`truncation_axis`.
Returns
-------
Refusal or None
The refusal, or ``None`` when the requested axis can be swept.
"""
from ..estimators.drtmle import DRTMLE
if refits_reduced_regressions(result) and isinstance(result.estimator, DRTMLE):
settings = result.estimator._reduction_configuration_refusal(
result.extra.get("drtmle"), result.nuisance
)
if settings is not None:
return replace(
settings,
condition=(
"a guarded DR-TMLE truncation curve refits the reduced regressions under "
"this estimator's reduction settings, and it refuses those settings: "
+ settings.condition
),
)
axis = truncation_axis(result, mechanism)
for _name, rule in _TRUNCATION_RULES:
refusal = rule(result, axis)
if refusal is not None:
return refusal
return None
[docs]
def truncation_curve(
result: TMLEResult,
bounds: Any = None,
*,
estimands: Any = None,
mechanism: bool | None = None,
) -> Any:
"""Re-estimate across a grid of truncation bounds.
Returns a tidy frame with one row per evaluated bound pair and estimand, giving the
point estimate and confidence interval. A fit whose inference status supplies no
inference, such as a selector-path collaborative fit, has no interval, so its frame
carries ``plugin_std_err``, ``plugin_interval_lower`` and
``plugin_interval_upper`` in place of ``std_err``, ``ci_lower`` and ``ci_upper``.
A scalar bound supplied by the caller remains
shorthand for the symmetric treatment-mechanism pair ``(bound, 1 - bound)``; an
observation- or intermediate-mechanism sweep uses ``(bound, 1)``. On an ordinary,
collaborative, or unguarded
doubly-robust fit only the targeting step is re-run -- the nuisance fits are cached
-- so each bound costs a small fraction of the original fit. A
``CTMLE(strategy="derivative")`` fit also re-runs its derivative update at each bound,
because both of its mechanisms are bounded there.
A **guarded DR-TMLE** fit is the exception, and it is why the combined report prices
this operation per result rather than per operation. Its targeting step alternates
against the reduced-dimension regressions, so every bound refits them. The ordinary
closure refits them at the fitted reduced bounds; the missing-outcome construction is
handed the swept bounds instead, because they define two of its regression targets.
The primary outcome regression and the propensity stay cached, so one bound costs
less than a fit -- but the default grid has several of them, so the **whole curve can
cost more than the fit it describes**. Budget it as refitting work.
A curve that is flat over the plausible range of bounds says the estimate does not
hinge on the truncation choice. A curve that drifts monotonically says the estimate
is sensitive to this finite-sample regularisation. It does not change the requested
estimand. The bound should be reported with the estimate.
Parameters
----------
result : TMLEResult
A fitted result.
bounds : sequence of float or None
Lower truncation values to try. ``None`` uses a grid from 0.001 to 0.2 and adds
every exact bound pair the requested parameters used. An explicit sequence is
evaluated exactly as requested and is never expanded with fitted points. On a
declared known treatment mechanism, a ``g(W)`` pair that moves a declared value is
one the fit refuses: an explicit one is refused, and the default grid omits it.
estimands : sequence of str or None
Restrict the emitted rows to these parameters. ``None`` uses every parameter the
fit reported. Two forms are accepted. A reported parameter selects its own row,
so ``"ey[high]"`` gives one row per bound pair. A registered target name selects
every reported alias of that target, in report order, so ``"ey"`` on a three-armed
fit gives three rows per bound pair and on a two-armed fit gives one. A name that
is neither is refused, because it has no fitted pair and no fitted estimate for a
row to reference, and an empty selection is refused because it names no row.
mechanism : bool or None
Sweep the bound on ``P(Delta = 1 | A, W)`` (and the intermediate density)
instead of the one on ``g(W)``. ``None`` resolves to the only axis the fit
has: the observation mechanism on a natural-course fit, which estimates no
treatment law, and ``g(W)`` on every other fit. That probability divides the
clever covariate exactly as the propensity does, so it has a truncation curve
for exactly the same reason -- and it is the one that goes unexamined, because
it has no familiar name. Requires a fit with ``delta=`` or ``intermediate=``.
An incremental fit refuses the ``g(W)`` axis, including the default, because its
estimand contains ``g``. ``truncation_refusal`` states every refusal, and
the capability row reads it too.
Note what the curve does and does not show. Truncating a mechanism cannot move
the *estimand*: the plug-in is an average of targeted predictions and contains
no mechanism at all. What moves is the second-order remainder, so a curve that
drifts is saying the estimate is leaning on rows the bound is holding up.
Returns
-------
dataframe
One row per evaluated bound pair and estimand, with the point estimate, interval,
fitted pair and estimate, and signed movement from that fitted estimate.
Raises
------
CapabilityError
When the requested axis is refused, as :func:`truncation_refusal` states, or when an
explicit ``g(W)`` bound moves a value of a declared known mechanism.
"""
estimator = result.estimator
if estimator is None:
raise CapabilityError(
"truncation_curve needs the fitted estimator that produced the result",
reason=RefusalReason.INVALID_DECLARATION,
)
# One predicate for the call and for the capability row, as ``tipping_gamma``'s
# ``use_ci=True`` has. The facade reads the same function to resolve its row, so a
# request the row calls available is a request this line admits.
refusal = truncation_refusal(result, mechanism)
if refusal is not None:
raise CapabilityError.from_refusal(refusal)
mechanism = truncation_axis(result, mechanism)
# Every reported alias against the target that answers it, in report order. The
# selection below reads it, and so does the row loop, so the registry is asked once.
reported_targets = {name: _target_name(result, name) for name in result.estimates}
reported = (
tuple(result.estimates) if estimands is None else _select(estimands, reported_targets)
)
def pair_for(lower: float) -> tuple[float, float]:
"""One lower bound as the pair the sweep evaluates, in the documented shorthand."""
return (lower, 1.0 if mechanism else 1.0 - lower)
# `retarget` takes *target* names, while `result.estimates` is keyed by the parameter
# names those targets reported -- the same thing for a two-armed fit, and `ey[high]`
# against `ey` for a wider one. The rows narrow back to `reported`: an arm-narrowed fit
# reports `ey[high]` alone while the estimator still retargets every arm of `ey`.
targets = {name: reported_targets[name] for name in reported}
names = tuple(dict.fromkeys(targets.values()))
pairs = {
name: _fitted_bound_pair(result, target, mechanism=mechanism)
for name, target in targets.items()
}
fitted_psi = {name: float(result.estimates[name].psi) for name in reported}
# A g(W) bound that moves a declared value is a bound the fit refuses, so the sweep
# never evaluates one. The default grid omits such a canned pair, and an explicit one
# is refused. The fitted pairs passed the fit's own check.
# The fit's own source decides: the declaration on the data, not the stored propensity,
# which a guarded DR-TMLE may replace with its fluctuated mechanism.
declared = None if mechanism else result.data.known_treatment
known = declared is not None
canned = (0.001, 0.005, 0.01, 0.025, 0.05, 0.1, 0.2)
if bounds is None:
admitted = (
pair
for pair in (pair_for(lower) for lower in canned)
if not known or _moved_known_count(result, pair) == 0
)
grid = sorted({*admitted, *pairs.values()})
else:
lowers = sorted(float(value) for value in bounds)
for lower in lowers:
if not 0.0 < lower < 0.5:
raise ValueError(f"truncation bounds must lie in (0, 0.5); got {lower}")
grid = [pair_for(lower) for lower in lowers]
for pair in grid if known else ():
count = _moved_known_count(result, pair)
if count and declared is not None:
refusal = moved_known_values_refusal(count, pair, declared.values)
raise CapabilityError.from_refusal(
replace(
refusal,
condition=(
"a truncation curve evaluates only the bound pairs that a fit "
"admits, and a fit refuses this one: " + refusal.condition
),
alternative=(
"Remove the bound from the grid, or pass bounds=None for the "
"default grid, which omits every pair that moves a known value."
),
)
)
rows: list[dict[str, Any]] = []
for lower, upper in grid:
if not 0.0 < lower < upper <= 1.0: # pragma: no cover - fitted config is validated
raise RuntimeError(f"invalid fitted truncation pair {(lower, upper)}")
pair = (lower, upper)
# Every draw, then combined the way the fit combines them. Sweeping one draw and
# calling the answer the fit's would compare a bound's effect on one split against
# a reported estimate that came from R -- and the difference between the two curves
# would read as sensitivity to the bound. Costs R times the sweep, which is still a
# fraction of one refit.
estimates = median_estimates(
[
estimator.retarget(
result.data,
repeat.nuisance,
estimands=names,
intermediate_value=result.intermediate_value,
g_bounds=None if mechanism else pair,
g_bounds_conditional=None if mechanism else pair,
nuisance_bound=lower if mechanism else None,
conditions=result.config.support_conditions,
)[0]
for repeat in result.repeats
],
)
# A retarget of `ey` answers every arm, and a fit narrowed to one arm asked for
# one row. Dropping the extra names is the ordinary case; missing a requested one
# is an internal inconsistency, so it is raised rather than silently omitted.
missing = [name for name in reported if name not in estimates]
if missing:
raise RuntimeError(
f"retargeting at bounds {pair} reported {sorted(estimates)}, which leaves "
f"{missing} of the requested parameters {list(reported)} unanswered"
)
# One truncation load per bound pair: it reads the pair and the mechanism flag,
# and nothing that varies between the parameters sharing that pair.
truncated_fraction = _clipped_fraction(result, pair, mechanism)
moved = _moved_factor(result, mechanism)
for name in reported:
estimate = estimates[name]
# The columns ``ParameterEstimate.to_dict`` publishes, for the same reason: a
# selector-path collaborative fit supplies no interval, so this frame reports
# its retained diagnostic under names that make no coverage claim rather than
# raising halfway through the sweep. One call names and computes all three,
# so the error and the interval cannot come from different branches.
fitted_lower, fitted_upper = pairs[name]
reference = fitted_psi[name]
rows.append(
{
"bound": lower,
"estimand": name,
"psi": estimate.psi,
**estimate.spread_columns(pvalue=False),
"truncated_fraction": truncated_fraction,
"is_fitted_bound": pair == (fitted_lower, fitted_upper),
# The upper bound and the fitted-reference columns follow the row's
# lower bound, estimate, spread and truncation columns.
"upper_bound": upper,
"fitted_lower_bound": fitted_lower,
"fitted_upper_bound": fitted_upper,
"fitted_psi": reference,
"delta_from_fitted": float(estimate.psi - reference),
# A composite fit divides by a product of factors, and a bound moves
# some of them only, so the row names the factors it truncated.
**({} if moved is None else {"truncated_factor": moved}),
**({"treatment_mechanism": "known"} if known else {}),
# The retarget classified the swept fit again, so an experimental
# configuration keeps its label and the sentence that names the condition.
**(
{"support": "experimental", "support_reason": estimate.support.reason}
if estimate.support.label == "experimental"
else {}
),
}
)
keys = list(dict.fromkeys(key for row in rows for key in row))
payload = {key: [row.get(key) for row in rows] for key in keys}
return result.data.frame_like(payload)
def _moved_known_count(result: TMLEResult, pair: tuple[float, float]) -> int:
"""How many units a ``g(W)`` bound pair moves on a declared known mechanism.
The mechanism is built from the declaration exactly as the fit's own check builds it,
and the rule is :meth:`Propensity.truncate
<cleverly.estimators._nuisance.Propensity.truncate>`, the one that check reads.
"""
known = result.data.known_treatment
if known is None: # pragma: no cover - the caller asks on a declared mechanism only
return 0
mechanism = KnownPropensity(np.asarray(known.values, dtype=float), result.data.arm_codes)
return int(np.count_nonzero(mechanism.truncate(pair).units))
def _select(estimands: Any, reported_targets: dict[str, str]) -> tuple[str, ...]:
"""The reported parameters an ``estimands=`` request selects, in report order.
A reported parameter selects its own row. A registered target name selects every
reported alias of that target, because the two spellings coincide only on a two-armed
fit: ``estimands=["ey"]`` has to keep working on a fit that reports ``ey[low]``,
``ey[medium]`` and ``ey[high]``, where the target is the name a caller has.
Every other name is refused here, before the first retarget. It has no fitted pair
and no fitted estimate, so every row it would fill is undefined. A typed arm label
the fit does not carry, such as ``ey[nope]``, is one of those: it is not stemmed back
to ``ey`` and answered with the whole target. An empty request is refused for the
same reason, and at the same point, because it names no row for the sweep to emit.
"""
requested = tuple(dict.fromkeys(estimands))
if not requested:
raise CapabilityError(
"estimands= selected no parameter; pass at least one of "
f"{list(reported_targets)}, or None for every parameter this fit reported",
reason=RefusalReason.INVALID_DECLARATION,
)
targets = [
name for name in dict.fromkeys(reported_targets.values()) if name not in reported_targets
]
unknown = [
name
for name in requested
if name not in reported_targets and name not in reported_targets.values()
]
if unknown:
also = f"; or a target name whose reported arms it selects: {targets}" if targets else ""
raise CapabilityError(
f"parameter(s) {unknown} were not reported by this fit; "
f"choose from {list(reported_targets)}{also}",
reason=RefusalReason.INVALID_DECLARATION,
)
selected = [
alias
for name in requested
for alias in (
(name,)
if name in reported_targets
else tuple(a for a, target in reported_targets.items() if target == name)
)
]
return tuple(dict.fromkeys(selected))
def _target_name(result: TMLEResult, name: str) -> str:
"""Return a reported alias's target, using its stem when no key was recorded."""
# A direct estimator result can leave the mapping empty. Its alias stem is then
# the only available target name.
keys = result.parameter_keys
key = keys.get(name) if keys else None
return key.estimand if key is not None else parameter_stem(name)
def _fitted_bound_pair(result: TMLEResult, target: str, *, mechanism: bool) -> tuple[float, float]:
"""The fitted mechanism pair applicable to one registered target."""
if mechanism:
return float(result.config.missingness_bound), 1.0
group = TARGETS[target].group
return g_bounds_for(group, result.config.g_bounds, result.config.g_bounds_conditional)
def _moved_factor(result: TMLEResult, mechanism: bool) -> str | None:
"""Which factors of a composite mechanism the swept bound truncates, or ``None``.
Only a composite-indicator fit divides by a product that a bound moves in part. The
``g_bounds`` axis truncates the treatment factor alone. The mechanism axis truncates
each observation factor. The composite itself is clipped inside the product of the
factor floors, which follows both.
"""
if result.extra.get("missing_data") != "composite":
return None
name = result.data.treatment_name
design = _design_factors(result)
if not mechanism:
if design:
return design[0]
conditioning = ", Delta_A=1" if result.data.has_missing_treatment else ""
return f"P({name}=a|W{conditioning})"
outcome_label, treatment_label = _observation_labels(result)
factors = []
if result.nuisance.missingness is not None:
factors.append(outcome_label if design else f"P(Delta=1|{name}=a,W)")
if result.nuisance.treatment_observation is not None:
factors.append(treatment_label)
return " and ".join(factors)
def _design_factors(result: TMLEResult) -> tuple[str, ...]:
"""The factor labels of a per-arm outcome-adaptive fit, or an empty tuple.
Such a fit regresses each factor on ``Qbar(a, W)``, not on ``W``, so the report names
the factor it fitted. :class:`~cleverly.estimators.ctmle.CTMLEOutcomeAdaptiveFit`
lists them, the treatment factor first.
"""
return tuple(getattr(result.ctmle_selection, "factors", ()))
def _composed_row(
result: TMLEResult, has_observation: bool, has_intermediate: bool, has_recording: bool
) -> str:
"""The name of a fit's derived denominator row.
A per-arm outcome-adaptive fit regresses each factor on ``Qbar(a, W)``, and inside each
baseline stratum with ``strata=``. Its product is the composite indicator's probability
given that design, so the row names the design's conditioning set, which the fit's
treatment factor label states. Every other fit names ``W``.
"""
name = _COMPOSED_ROWS[has_observation, has_intermediate, has_recording]
design = _design_factors(result)
if not design or has_intermediate:
return name
condition = _DESIGN_CONDITIONS[1] if design[0].endswith(",S)") else _DESIGN_CONDITIONS[0]
return _DESIGN_ROWS[name, condition]
def _observation_labels(result: TMLEResult) -> tuple[str, str]:
"""The labels of the outcome and of the treatment observation factor of a fit."""
design = _design_factors(result)
outcome = next((f for f in design if f.startswith("P(Delta=1|")), "P(Delta=1|A,W)")
treatment = next((f for f in design if f.startswith("P(Delta_A=1|")), TREATMENT_OBSERVATION)
return outcome, treatment
def _clipped_fraction(result: TMLEResult, pair: tuple[float, float], mechanism: bool) -> float:
"""Share of the nuisance the bound would clip, for whichever bound is swept.
The treatment branch asks :meth:`~cleverly.estimators._nuisance.Propensity.truncate`
and reports the share of **units** it moves at ``pair``. A unit is one row of the
mechanism, and one binding denominator is enough to make that row use a bounded
rather than fitted mechanism value. The share therefore belongs to the pair it was
evaluated at, and it
equals the :func:`positivity_report` figure only at that report's own pair. A
continuous treatment fits no arms at all, and the share is ``nan`` there rather than
a claim that the bound moved nothing.
The mechanism branch has one column per arm as well, but its bound is applied cell by
cell and is reported that way.
"""
if not mechanism:
return result.nuisance.propensity.truncate(pair).fraction
lower, upper = pair
# The intermediate entry must be the density for the level being targeted, not the
# raw P(Z = 1 | A, W): at z = 0 the covariate divides by the complement, so reading
# the array directly counts the wrong tail and reports a mirror-inverted fraction.
# This column is not a nicety -- it is the only part of the mechanism truncation
# curve that can detect a q_z positivity problem at all. A density that is clipped
# to a constant rescales both clever-covariate columns by the same factor, and the
# fluctuation ``epsilon * h`` is invariant to that, so ``psi`` sits flat across the
# whole sweep however badly the bound is binding.
missingness = (
result.nuisance.missingness_at_realised_arm(result.data.treatment)
if not result.nuisance.fits_treatment
else result.nuisance.missingness
)
candidates = [missingness, result.nuisance.treatment_observation]
if result.nuisance.intermediate is not None and result.intermediate_value is not None:
candidates.append(result.nuisance.intermediate_density(result.intermediate_value, 0.0))
parts = [
np.asarray(values, dtype=float).reshape(-1) for values in candidates if values is not None
]
divisor = np.concatenate(parts)
return float(np.mean((divisor < lower) | (divisor > upper)))