r"""E-values (VanderWeele & Ding, 2017).
The E-value answers the same question as the omitted-variable bound but on the risk
ratio scale, and with a different parameterisation that has become the convention in
epidemiology:
*the minimum strength, on the risk ratio scale, that each of two confounder
associations would need if both had equal strength: one with treatment and one
with outcome, conditional on the measured covariates, to explain away the
observed association.*
Unequal associations can trade off: one can be below the E-value if the other is
above it. The E-value alone does not establish whether such confounding is plausible.
For an observed risk ratio :math:`RR \ge 1`,
.. math:: E = RR + \sqrt{RR\,(RR - 1)}
and for :math:`RR < 1` the same formula is applied to :math:`1/RR`. An E-value of 1
means no unmeasured confounding is needed -- the estimate is already
compatible with the null. A larger E-value raises the equal-strength threshold.
The E-value for the *confidence limit* is usually the more important number: it says
how strong a confounder would need to be to move the interval to include the null,
which is the claim a sceptical reader will make.
Conversions for other effect scales are approximations and are flagged as such:
* **odds ratio**, common outcome: :math:`RR \approx \sqrt{OR}` (VanderWeele, 2017).
This approximation can lie above or below the risk ratio. For rare outcomes,
the odds ratio itself approximates the risk ratio.
* **continuous outcome**: standardise the difference by the outcome's standard
deviation. Chinn's :math:`d \approx \log(OR) / 1.81` gives an odds ratio; the
additional common-outcome conversion :math:`RR \approx \sqrt{OR}` then gives
:math:`RR \approx \exp(0.905 \times d)`. This assumes approximate normality and is
only a rough guide -- prefer
:func:`cleverly.sensitivity.omitted_variable_bounds` for a continuous outcome,
which needs no such conversion.
This package has no implemented E-value bound for its fitted controlled direct effect
and confounding model. :func:`_select_evalue` refuses every request on a fit with an
intermediate variable, because the ordinary conversion does not establish sensitivity
to unmeasured confounding for that target.
"""
from __future__ import annotations
from dataclasses import dataclass, replace
from typing import TYPE_CHECKING, Any, Literal
import numpy as np
from .._inference_status import status_record
from .._refusal import NO_ALTERNATIVE, Refusal, RefusalReason
from .._support import UNCLASSIFIED, Support, combine_support, record_refusal_site
from ..assessment import AssessmentStatus
from ..estimators.direct_effect import declares_intermediate
from ..exceptions import CapabilityError
from ..interventions.learned import LEARNED_RULE_EVALUE_REFUSAL
from ..targets.population_intervention import (
NATURAL_COURSE_TARGET,
POPULATION_INTERVENTION_TARGETS,
)
from ._derived import _derived_risk_ratio, _risk_ratio_refusal
from ._parameters import arm_parameter_keys
if TYPE_CHECKING: # pragma: no cover - typing only
from ..estimators.base import TMLEResult
from ..inference.influence import ParameterEstimate
__all__ = ["EValue", "evalue", "evalue_from_rr"]
#: Chinn's OR-to-SMD factor followed by the common-outcome square-root conversion.
_SMD_TO_LOG_RR = 1.81 / 2
#: Why the standardized conversion stops on a fit with a response mechanism. It names
#: the quantity the conversion lacks, which is the *population* standard deviation of the
#: outcome, and it says which E-value branches stay open.
_STANDARDISED_MISSING_REFUSAL = Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"the standardized E-value has no population outcome standard deviation for a fit "
"with a response mechanism. The conversion divides the difference by sd(Y) computed "
"from the observed rows alone, and under missing at random the respondents' standard "
"deviation estimates a different quantity from the population standard deviation the "
"estimate is standardized against.",
alternative=(
"An E-value computed from a reported risk ratio or odds ratio reads only the "
"estimate and its interval, so it stays available."
),
roadmap="X65",
)
#: Why every E-value branch stops on a fit that supplies no inference. Each branch reads
#: the estimate's interval, or the reference arm's standard error, and such a fit supplies
#: neither. Raised in :func:`_select_evalue` rather than where the interval is read, so
#: the capability row reports ``unavailable`` instead of the computation raising under a
#: row that advertised ``available=True``. The status's own reason follows it as a new
#: sentence.
_EVALUE_NEEDS_INFERENCE = "an E-value is built from the reported estimate and its interval. "
#: Why every E-value branch stops on a fit with an intermediate variable. This package
#: has no implemented bound for its fitted controlled direct effect and confounding model.
#: The E-value section of ``docs/technical-reference/validation-methods.md`` gives the
#: support boundary, and X49 in ``docs/roadmap.md`` holds the per-arm bound that fills the
#: gap. Raised in
#: :func:`_select_evalue` before the
#: estimand is resolved, so the row reads ``unavailable`` on every such fit with a discrete
#: treatment, a multi-arm fit included, and no branch computes. A continuous-treatment fit
#: meets the check before it and reads ``not_applicable``, so it is refused as well.
_DIRECT_EFFECT_REFUSAL = Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"this package has no implemented E-value bound for the fitted controlled direct "
"effect and its confounding model. Applying the ordinary E-value conversion does "
"not establish sensitivity to unmeasured treatment-outcome or intermediate-outcome "
"confounding for this target, because its bias splits into one bounding factor per "
"arm and the ordinary conversion inverts one factor. docs/roadmap.md X49 tracks the "
"E-value of their product.",
alternative=NO_ALTERNATIVE,
roadmap="X49",
)
def _evalue_inference_refusal(status: str) -> Refusal:
"""The E-value refusal at a non-inferential status, built on that status's reason."""
record = status_record(status)
record_refusal_site("evalue.inference_status", RefusalReason.INFERENCE_GAP, record.reopened_by)
return Refusal(
RefusalReason.INFERENCE_GAP,
_EVALUE_NEEDS_INFERENCE + record.reason,
alternative=record.alternative,
roadmap=record.reopened_by,
)
def _refuse(status: AssessmentStatus, refusal: Refusal) -> _EValueRefusal:
"""Build the E-value refusal that carries one complete record."""
return _EValueRefusal(
status,
refusal.condition,
reason=refusal.reason,
alternative=refusal.alternative,
roadmap=refusal.roadmap,
)
def evalue_from_rr(risk_ratio: float) -> float:
"""E-value for a risk ratio.
>>> round(evalue_from_rr(2.0), 4)
3.4142
>>> evalue_from_rr(1.0)
1.0
>>> round(evalue_from_rr(0.5), 4) # symmetric under inversion
3.4142
"""
rr = float(risk_ratio)
if not np.isfinite(rr) or rr <= 0:
return float("nan")
if rr < 1.0:
rr = 1.0 / rr
return float(rr + np.sqrt(rr * (rr - 1.0)))
def _evalue_for_limit(limit: float, *, above_null: bool) -> float:
"""E-value for a confidence limit, which is 1 when the limit crosses the null.
The null test runs before any positivity test. A limit outside the risk-ratio
parameter space still says which side of the null the interval reaches, and an
interval that already covers the null needs no confounding to explain it. The
fixed-baseline conversion is affine, so a risk-difference limit below the negative of
the baseline risk maps to a nonpositive ratio limit. A point ratio above the null
reads that lower limit, and a nonpositive value there reports 1, not ``nan``. Only
``nan`` itself, which names no side of the null, reports ``nan``.
A point ratio below the null reads the *upper* limit instead, and the fixed-baseline
conversion never sends that bound out of the parameter space:
:func:`~cleverly.inference.delta.wald_ci` gives ``high >= psi``, and
:func:`_reject_unusable_baseline` has already refused a fit whose ``baseline.psi`` is
nonpositive or whose ``baseline.psi + psi`` is nonpositive.
"""
if np.isnan(limit):
return float("nan")
if above_null:
return 1.0 if limit <= 1.0 else evalue_from_rr(limit)
return 1.0 if limit >= 1.0 else evalue_from_rr(limit)
[docs]
@dataclass(frozen=True)
class EValue:
"""E-values for a point estimate and the confidence limit nearest the null.
Parameters
----------
estimand : str
Alias of the estimand this row describes.
scale : str
Scale the estimate was reported on before conversion.
risk_ratio : float
The estimate expressed as a risk ratio.
risk_ratio_ci : tuple of float
That ratio's confidence interval. A converted lower bound below zero is outside
the risk-ratio parameter space. The report truncates it at zero, records the
untruncated value in :attr:`truncated_bound`, and says so in :attr:`note`.
point : float
Equal-strength threshold for the confounder's treatment and outcome
associations on the risk-ratio scale to explain the point estimate away.
Unequal association strengths can trade off.
limit : float
The equal-strength threshold to move the confidence limit to the null.
approximate : bool
Whether reaching the risk-ratio scale needed an approximate conversion.
note : str
Explanation of that conversion, when one was used.
source_estimand : str or None
Reported source alias, which can differ from a derived output alias.
truncated_bound : float or None
The converted lower bound before truncation, when the report truncated one.
``None`` when :attr:`risk_ratio_ci` reports the conversion unchanged.
support : Support
The support of the source estimate, combined with the derived ratio's when the
row derives one. An E-value built outside a classified fit is experimental.
"""
estimand: str
scale: str
risk_ratio: float
risk_ratio_ci: tuple[float, float]
point: float
limit: float
approximate: bool
note: str = ""
source_estimand: str | None = None
truncated_bound: float | None = None
support: Support = UNCLASSIFIED
[docs]
def summary(self) -> str:
"""Return a printable summary.
Returns
-------
str
A printable report, one line per reported quantity.
"""
marker = " [experimental]" if self.support.label == "experimental" else ""
lines = [
f"E-value for {self.estimand!r} ({self.scale} scale){marker}",
"-" * 40,
f"risk ratio {self.risk_ratio:.4g} "
f"[{self.risk_ratio_ci[0]:.4g}, {self.risk_ratio_ci[1]:.4g}]"
+ (" (approximate conversion)" if self.approximate else ""),
f"E-value, point estimate : {self.point:.4f}",
f"E-value, confidence limit : {self.limit:.4f}",
"",
f"If its two risk-ratio associations had equal strength, an unmeasured "
f"confounder would need {self.point:.2f} with both treatment and outcome, "
f"beyond the measured covariates, to explain away the point estimate; "
+ (
f"{self.limit:.2f} to move the interval to the null."
if self.limit > 1.0
else "the interval already includes the null."
),
"Unequal strengths can trade off: one association may be below the E-value "
"if the other is above it.",
]
if self.note:
lines.extend(["", self.note])
if self.support.label == "experimental":
lines.extend(["", f"experimental: {self.support.reason}"])
return "\n".join(lines)
[docs]
def to_dict(self) -> dict[str, float | str | bool | None]:
"""Return a JSON-compatible representation.
``truncated_bound`` carries the untruncated lower bound, so a caller who reads
only this mapping still sees that ``rr_ci_lower`` is a boundary value rather than
a converted confidence limit.
Returns
-------
dict
A JSON-compatible mapping of every reported field.
"""
return {
"estimand": self.estimand,
"scale": self.scale,
"risk_ratio": self.risk_ratio,
"rr_ci_lower": self.risk_ratio_ci[0],
"rr_ci_upper": self.risk_ratio_ci[1],
"evalue_point": self.point,
"evalue_limit": self.limit,
"approximate": self.approximate,
"note": self.note,
"source_estimand": self.source_estimand,
"truncated_bound": self.truncated_bound,
"support": self.support.label,
}
def __repr__(self) -> str: # pragma: no cover - debugging aid
return self.summary()
class _EValueRefusal(CapabilityError):
"""One refusal, carrying the status :mod:`cleverly.assessment` publishes for it.
The status is the enum member rather than its string, so every branch below is checked
where it is written. Thirteen call sites spelled the value as a bare literal and only
the one place that built an :class:`~cleverly.assessment.AssessmentStatus` from it
validated any of them, which made a misspelling a runtime error on the row it refused.
The import is module level and introduces no cycle. :mod:`cleverly.assessment` reaches
this module through function-local imports only, and nothing on its own module-level
import chain reads :mod:`cleverly.sensitivity`.
"""
def __init__(
self,
status: AssessmentStatus,
text: str,
/,
*,
reason: RefusalReason,
alternative: str | None = None,
roadmap: str | None = None,
) -> None:
super().__init__(text, reason=reason, alternative=alternative, roadmap=roadmap)
self.status = status
#: The conversions this module implements. Each one has a branch in
#: :func:`_evalue_from_selection`, and :mod:`cleverly.assessment` reads two of these
#: strings to decide what a combined run costs.
_Branch = Literal[
"reported_rr",
"reported_or",
"derived_rr",
"fixed_baseline_ate",
"gaussian_difference",
]
#: The one branch that has to retarget cached nuisances rather than read a stored
#: estimate. :mod:`cleverly.assessment` prices a combined run by it and decides by it
#: whether an argument-free request replays under its own estimand name, so the string is
#: named here rather than repeated there.
_DERIVED_RR: _Branch = "derived_rr"
@dataclass(frozen=True)
class _EValueSelection:
source: str
branch: _Branch
def _default_estimand(result: TMLEResult, keys: dict[str, Any]) -> str:
"""Choose a reported contrast from structured parameter identities."""
for target in ("rr", "or", "ate", "att", "atc"):
candidates = [
name
for name in result.estimates
if (key := keys.get(name)) is not None
and key.estimand == target
and key.axis == "arm"
and key.reference is not None
and key.stratum is None
]
if len(candidates) > 1:
raise _EValueRefusal(
AssessmentStatus.DEFERRED,
f"an E-value needs one contrast; choose an explicit estimand from {candidates}",
reason=RefusalReason.INVALID_DECLARATION,
)
if candidates:
return candidates[0]
raise _EValueRefusal(
AssessmentStatus.NOT_APPLICABLE,
"an E-value needs a supported two-arm contrast; this fit reports no such parameter",
reason=RefusalReason.INVALID_DECLARATION,
)
def _baseline_mean(result: TMLEResult, estimand: str, keys: dict[str, Any]) -> str | None:
"""The reported mean under the arm ``estimand`` is contrasted against, if any.
The risk-difference conversion divides by the baseline risk, and the baseline is the
contrast's *reference* arm rather than arm ``0``: with ``reference=1`` a two-armed
fit's ``ate`` is ``E[Y^0] - E[Y^1]`` and dividing by ``EY0`` would convert a
difference into a ratio of the wrong pair. Found by arm through
:func:`~cleverly.sensitivity._parameters.arm_parameter_keys`, so ``ey0`` and
``ey[low]`` are both recognised -- whichever of them the fit was asked for.
"""
parameter = keys.get(estimand)
if parameter is None or parameter.reference is None:
return None
return next(
(
name
for name, candidate in keys.items()
if name in result.estimates
and candidate.axis == "arm"
and candidate.estimand in {"ey", "ey0", "ey1"}
and candidate.reference is None
and candidate.value == parameter.reference
and candidate.stratum is None
),
None,
)
def _select_evalue(result: TMLEResult, estimand: str | None) -> _EValueSelection:
"""Select the exact, approximate, or refusal branch for one request."""
if getattr(result, "assessment_family", None) != "point":
raise _EValueRefusal(
AssessmentStatus.UNAVAILABLE,
"no longitudinal sensitivity derivation is registered for an E-value: the "
"bounding factor needs a per-period latent law for a sequential g-formula "
"target (docs/roadmap.md F16 tracks this stop).",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=NO_ALTERNATIVE,
roadmap="F16",
)
if not result.data.has_arms:
raise _EValueRefusal(
AssessmentStatus.NOT_APPLICABLE,
"an E-value requires a discrete arm contrast",
reason=RefusalReason.INVALID_DECLARATION,
)
if result.config.parameter_axis == "learned_rule":
raise _EValueRefusal(
AssessmentStatus.UNAVAILABLE,
"an E-value does not cover a learned-rule fit: " + LEARNED_RULE_EVALUE_REFUSAL + ".",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=NO_ALTERNATIVE,
roadmap="F27",
)
# Fit-wide and ahead of the estimand, so a level or a multi-arm default refuses here too.
if declares_intermediate(result):
raise _refuse(AssessmentStatus.UNAVAILABLE, _DIRECT_EFFECT_REFUSAL)
keys = arm_parameter_keys(result)
explicit = estimand is not None
source = _default_estimand(result, keys) if estimand is None else estimand
if source not in result.estimates:
raise _EValueRefusal(
AssessmentStatus.UNAVAILABLE,
f"estimand {source!r} was not requested in this fit",
reason=RefusalReason.INVALID_DECLARATION,
)
key = keys.get(source)
attributable = source.split("[", 1)[0] if key is None else key.estimand
if attributable in POPULATION_INTERVENTION_TARGETS:
raise _EValueRefusal(
AssessmentStatus.NOT_APPLICABLE,
"an E-value is defined for a two-arm risk ratio or contrast, and "
f"{attributable} compares the natural course with one reference arm.",
reason=RefusalReason.DIFFERENT_ESTIMAND,
alternative=(
"Call evalue(result, 'rr') or evalue(result, 'ate') on a fit that also "
"requests the two-arm contrast, such as estimands=('par', 'ate', 'rr')."
),
)
if attributable == NATURAL_COURSE_TARGET:
raise _EValueRefusal(
AssessmentStatus.NOT_APPLICABLE,
"an E-value is defined for a two-arm risk ratio or contrast, and ey_obs is the "
"natural-course mean, not a two-arm contrast.",
reason=RefusalReason.DIFFERENT_ESTIMAND,
alternative=(
"Call evalue(result, 'rr') or evalue(result, 'ate') on a fit that also "
"requests the two-arm contrast."
),
)
if key is None:
raise _EValueRefusal(
AssessmentStatus.UNAVAILABLE,
f"estimand {source!r} has no structured parameter key",
reason=RefusalReason.INVALID_DECLARATION,
)
if key.axis != "arm":
raise _EValueRefusal(
AssessmentStatus.NOT_APPLICABLE,
f"an E-value needs an arm contrast, not the {key.axis!r} axis",
reason=RefusalReason.INVALID_DECLARATION,
)
if key.reference is None:
raise _EValueRefusal(
AssessmentStatus.NOT_APPLICABLE,
f"an E-value needs a two-arm contrast; target {key.estimand!r} has no reference arm",
reason=RefusalReason.INVALID_DECLARATION,
)
if key.stratum is not None:
raise _EValueRefusal(
AssessmentStatus.NOT_APPLICABLE,
"this package requires an unconditioned marginal arm contrast; this contrast "
"is conditional on a baseline stratum.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=(
"Call evalue on the marginal alias of the same fit, or fit the stratum's "
"rows as their own study."
),
roadmap="X65",
)
# Fit-wide rather than branch-scoped, unlike the standardized rule below: the reported
# and derived ratio branches read ``estimate.ci`` and the Gaussian branch reads the
# reference arm's ``std_error``, so no branch survives a fit that supplies neither.
# It runs after the request checks directly above, the requested name, its key and its
# axis, which say that an E-value is not defined for this request on any fit. So on a
# fit without an intermediate variable, a level such as ``ey1`` stays ``not_applicable``
# here too.
if not result.estimates[source].supplies_inference:
raise _refuse(
AssessmentStatus.UNAVAILABLE,
_evalue_inference_refusal(result.estimates[source].inference),
)
if key.estimand == "rr":
return _EValueSelection(source, "reported_rr")
if key.estimand == "or" and explicit:
return _EValueSelection(source, "reported_or")
if result.config.family == "gaussian" and key.estimand in {"ate", "att", "atc"}:
# Raised here, in the selection, rather than in the execution below. The
# capability row is built from this function, so a refusal placed in
# ``_evalue_from_selection`` would leave the row advertising ``available=True``
# beside a call that raises. Branch-scoped rather than fit-wide: the reported and
# derived ratio branches read no standard deviation at all.
if result.data.has_missing_outcome:
raise _refuse(AssessmentStatus.UNAVAILABLE, _STANDARDISED_MISSING_REFUSAL)
return _EValueSelection(source, "gaussian_difference")
if key.estimand in {"att", "atc"}:
raise _EValueRefusal(
AssessmentStatus.UNAVAILABLE,
f"{key.estimand.upper()} needs a conditional baseline risk and a supported "
"conditional ratio target; neither is available.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative=(
"Call evalue(result, 'ate') for the marginal ATE of the same arms, which "
"describes a different population. On a gaussian outcome the standardized "
"ATT or ATC E-value runs."
),
roadmap="F34",
)
if key.estimand not in {"ate", "or"}:
raise _EValueRefusal(
AssessmentStatus.NOT_APPLICABLE,
f"an E-value has no supported two-arm contrast for target {key.estimand!r}",
reason=RefusalReason.INVALID_DECLARATION,
)
refusal = _risk_ratio_refusal(result, source, keys)
if refusal is None:
return _EValueSelection(source, "derived_rr")
if key.estimand == "or":
return _EValueSelection(source, "reported_or")
baseline = _baseline_mean(result, source, keys)
if baseline is None or not np.isfinite(result[baseline].psi) or result[baseline].psi <= 0:
raise _refuse(
AssessmentStatus.UNAVAILABLE,
replace(
refusal,
condition=refusal.condition.rstrip(".")
+ "; a finite positive reported reference-arm mean is also absent.",
),
)
_reject_unusable_baseline(result[baseline], result[source], refusal)
return _EValueSelection(source, "fixed_baseline_ate")
def _reject_unusable_baseline(
baseline: ParameterEstimate, estimate: ParameterEstimate, refusal: Refusal
) -> None:
"""Refuse a fixed-baseline conversion whose ratio leaves the parameter space.
The conversion holds the reference risk fixed and divides by it, so a positive sign
is not enough. A reference risk its own standard error does not separate from zero
is an unstable denominator: it sends the ratio to any value the caller cares to name.
A risk difference at or below the negative of that risk implies a nonpositive risk in
the contrast arm, which is outside the risk-ratio parameter space. Both are refused
by name rather than reported as an ordinary approximate E-value.
"""
# The inner record keeps its code, alternative and roadmap id: the fallback adds a
# reason why the approximate route fails too, and the exact route stays the remedy.
error = baseline.std_error
if not np.isfinite(error) or baseline.psi <= error:
raise _refuse(
AssessmentStatus.UNAVAILABLE,
replace(
refusal,
condition=f"{refusal.condition.rstrip('.')}; the reported reference-arm mean "
f"{baseline.psi:.4g} is not separated from zero by its own standard error "
f"{error:.4g}, so the fixed-baseline conversion has no stable denominator.",
),
)
ratio = (baseline.psi + estimate.psi) / baseline.psi
if not np.isfinite(ratio) or ratio <= 0:
raise _refuse(
AssessmentStatus.UNAVAILABLE,
replace(
refusal,
condition=f"{refusal.condition.rstrip('.')}; the reported difference "
f"{estimate.psi:.4g} against the reference-arm mean {baseline.psi:.4g} "
"implies a nonpositive risk in the contrast arm, which no risk ratio "
"describes.",
),
)
[docs]
def evalue(result: TMLEResult, estimand: str | None = None) -> EValue:
"""E-value for a fitted result.
``estimand=None`` prefers a reported risk ratio. For an eligible binomial
ordinary-TMLE result that reports only a marginal ATE or odds ratio, it retargets the
retained nuisances to the matching exact risk ratio. An explicit odds-ratio request
preserves the caller's chosen square-root approximation. Routing uses structured
parameter keys or forward-composed fitted arm identities; display aliases are never parsed.
Unsupported exact retargets can use a reported reference risk for an approximate ATE
conversion. That conversion needs a reference risk its own standard error separates
from zero, and it refuses a difference that implies a nonpositive risk in the contrast
arm. Gaussian ATE, ATT, and ATC use the standardized-difference approximation, which
divides by the weighted outcome standard deviation on a weighted fit.
Parameters
----------
result : TMLEResult
A fitted result.
estimand : str or None
Alias to report on. ``None`` picks the ratio-scale estimand the fit
reported, which is the scale an E-value is defined on.
Returns
-------
EValue
The association needed to explain the estimate and the interval away.
Raises
------
CapabilityError
If no branch covers the request. A fit with an intermediate variable estimates a
controlled direct effect, and every request on it refuses because this package
has no implemented E-value bound for its fitted target and confounding model.
"""
return _evalue_from_selection(result, _select_evalue(result, estimand))
def _standardising_sd(result: TMLEResult) -> tuple[float, bool]:
"""The outcome standard deviation the Gaussian conversion divides by.
A weighted fit estimates a weighted parameter, so the standard deviation that
standardises it belongs to the population the observation weights describe. The
unweighted sample deviation would put the numerator and the denominator on two
scales. The weighted form carries the reliability-weight correction
``sum(w) - sum(w^2) / sum(w)``, which is ``n - 1`` when every weight is one. So an
unweighted fit reports the plain sample deviation.
Parameters
----------
result : TMLEResult
A fitted result whose outcome is observed on at least two rows.
Returns
-------
tuple of (float, bool)
The standard deviation, and whether the observation weights moved it.
Raises
------
CapabilityError
If the correction is not positive, which names a row count rather than a
variance. One observed row gives ``w - w^2 / w = 0``, and no observed row of
positive weight gives ``0`` as well.
"""
observed = result.data.outcome[result.data.observed]
weights = np.asarray(result.data.weights, dtype=float)[result.data.observed]
total = float(np.sum(weights))
correction = total - float(np.sum(weights**2)) / total if total > 0 else 0.0
if correction <= 0:
raise CapabilityError(
"cannot standardise the effect: a standard deviation needs at least two "
f"observed rows of positive weight, and the reliability correction is "
f"{correction:.4g}",
reason=RefusalReason.INVALID_DECLARATION,
)
mean = float(np.average(observed, weights=weights))
variance = float(np.sum(weights * (observed - mean) ** 2)) / correction
return float(np.sqrt(variance)), bool(result.data.is_weighted)
def _evalue_from_selection(result: TMLEResult, selection: _EValueSelection) -> EValue:
"""Execute the same branch used to declare this request available."""
source = selection.source
estimate = result[source]
low, high = estimate.ci
reported_estimand = source
truncated_bound: float | None = None
support = estimate.support
if selection.branch == "reported_rr":
rr, ci, approximate, note, scale = (
estimate.psi,
(low, high),
False,
"",
"risk ratio",
)
elif selection.branch == "reported_or":
rr = float(np.sqrt(estimate.psi))
ci = (float(np.sqrt(low)), float(np.sqrt(high)))
approximate = True
scale = "odds ratio"
note = (
"Converted with RR ~ sqrt(OR), an approximation for common outcomes. "
"It can lie above or below the risk ratio; prefer a directly estimated risk ratio."
)
elif selection.branch == "derived_rr":
derived = _derived_risk_ratio(result, source)
rr, ci, approximate, scale = derived.psi, derived.ci, False, "risk ratio"
reported_estimand = derived.name
support = combine_support([estimate.support, derived.support], kind="derived")
note = (
f"Derived {derived.name!r} from cached nuisances by retargeting the reported "
f"source contrast {source!r}; no nuisance model was refit."
)
elif selection.branch == "fixed_baseline_ate":
baseline_name = _baseline_mean(result, source, arm_parameter_keys(result))
assert baseline_name is not None
baseline = result[baseline_name]
rr = float((baseline.psi + estimate.psi) / baseline.psi)
raw = (
float((baseline.psi + low) / baseline.psi),
float((baseline.psi + high) / baseline.psi),
)
# Only the lower bound can leave the parameter space. ``wald_ci`` gives
# ``high >= psi`` on a difference scale, and ``_reject_unusable_baseline`` has
# already refused this fit unless ``baseline.psi > 0`` and
# ``baseline.psi + estimate.psi > 0``, so ``raw[1] >= rr > 0``.
truncated_bound = raw[0] if raw[0] < 0.0 else None
ci = (0.0, raw[1]) if truncated_bound is not None else raw
approximate, scale = True, "risk difference"
note = (
"Converted from a risk difference using the reported reference-arm risk. "
"The conversion holds that risk fixed and ignores its sampling error."
)
if truncated_bound is not None:
note += (
" The difference interval reaches below the negative of that risk, so the "
f"converted lower bound {truncated_bound:.4g} lies outside the risk-ratio "
"parameter space. The report truncates it at 0, the smallest risk ratio "
"there is, and keeps the untruncated value in 'truncated_bound'. "
)
note += (
"The point ratio is above the null and the lower bound is below it, so "
"the interval covers the null and the confidence-limit E-value is 1."
if rr >= 1.0
else "The point ratio is below the null, so the confidence-limit E-value "
"reads the upper bound, which the truncation does not touch. That "
"interval can still exclude the null."
)
elif selection.branch == "gaussian_difference":
sd, weighted = _standardising_sd(result)
if not np.isfinite(sd) or sd <= 0:
raise CapabilityError(
"cannot standardise the effect: the outcome has no variance",
reason=RefusalReason.INVALID_DECLARATION,
)
rr = float(np.exp(_SMD_TO_LOG_RR * estimate.psi / sd))
ci = (
float(np.exp(_SMD_TO_LOG_RR * low / sd)),
float(np.exp(_SMD_TO_LOG_RR * high / sd)),
)
approximate, scale = True, "mean difference"
note = (
f"Standardised by {'weighted ' if weighted else ''}sd(Y) = {sd:.4g}. "
"Chinn's log(OR) / 1.81 step maps the "
"standardised mean difference to an odds ratio. The subsequent common-outcome "
"square-root conversion and E-value calculation are additional approximations."
)
else:
raise AssertionError(selection.branch)
if selection.branch in ("reported_or", "fixed_baseline_ate", "gaussian_difference"):
# A converted effect is a derived estimate, so the row must list the derived kind.
support = combine_support([support], kind="derived")
above_null = rr >= 1.0
limit = ci[0] if above_null else ci[1]
return EValue(
estimand=reported_estimand,
scale=scale,
risk_ratio=rr,
risk_ratio_ci=ci,
point=evalue_from_rr(rr),
limit=_evalue_for_limit(limit, above_null=above_null),
approximate=approximate,
note=note,
source_estimand=source,
truncated_bound=truncated_bound,
support=support,
)