r"""Estimands and the influence curves their estimates are reported from.
Efficient influence curves, everywhere but one: :func:`reduced_corrections` builds the two
terms a doubly-robust fit *subtracts* from :math:`D^*`, and the difference is the
estimator's own asymptotic influence function at the nuisance limits rather than the
canonical gradient at :math:`P_0`. The distinction is set out there and in
:mod:`cleverly.estimators.drtmle`; it is called out in this first line because a reader who
knows that TMLE's curve is the efficient one has no other reason to think this module holds
anything else.
Most single-draw estimates report a weighted targeted mean and the centered sample variance of
their influence curve divided by :math:`n`. Two exceptions exist. The stacked missing-outcome
natural-course CV-TMLE uses the curve's raw second moment. Fold-evaluated CV-TMLE uses the
cross-validated variance of the fold curves. Repeated cross-fitting reports the median point and
median split-adjusted variance while retaining a central draw's curve for marginal diagnostics.
Writing the curve explicitly supports cluster-robust variance, simultaneous bands, the delta
method, and score diagnostics where the declared covariance rule defines them.
The influence curves, on the ``[0, 1]`` outcome scale, with
:math:`r_i = \Delta_i (Y_i - \bar Q^*(A_i, W_i))`:
.. math::
\mathrm{IC}^{EY_1}_i &= h_1(A_i, W_i)\, r_i + \bar Q^*(1, W_i) - \psi_1 \\
\mathrm{IC}^{EY_0}_i &= h_0(A_i, W_i)\, r_i + \bar Q^*(0, W_i) - \psi_0 \\
\mathrm{IC}^{ATE}_i &= \mathrm{IC}^{EY_1}_i - \mathrm{IC}^{EY_0}_i \\
\mathrm{IC}^{ATT[a]}_i &= h_a(A_i, W_i)\, r_i
+ \frac{\mathbb 1\{A_i = a\}}{P(A = a)}
\bigl(\bar Q^*(a, W_i) - \bar Q^*(r, W_i) - \psi_{a}\bigr)
with :math:`r` the reference arm, one such curve per non-reference arm, and the ATC
mirrors it with :math:`\mathbb 1\{A_i = r\} / P(A = r)` in place of the arm's own
indicator -- every one of its parameters conditioning on the same population. With two
arms this is the classic ATT, :math:`a = 1` and :math:`r = 0`. Note that these curves
carry an extra term beyond "clever covariate times residual": the estimand conditions on a
*random* event (``A = a``), so the uncertainty in ``P(A = a)`` contributes. Omitting
it -- a common bug -- understates the standard error.
Observation weights enter as :math:`\mathrm{IC}_i \mapsto \tilde w_i \mathrm{IC}_i` with
:math:`\tilde w` normalised to mean one, so weighted and unweighted variances are
directly comparable. That row-wise multiplication is not a convenience: with every
nuisance fitted and every average taken under the weighted empirical measure, the
estimand is the same functional evaluated at the tilted law
:math:`dP_w = w\,dP / E[w]`, and its efficient influence function is exactly
:math:`(w / E[w])\, D^*_{P_w}`. :mod:`cleverly.data.weighting` derives that, states
which weighting problems it does and does not cover, and
``tests/unit/test_weighted_estimand.py`` verifies it numerically against a longhand
statement of :math:`\Psi(P_w)`.
"""
from __future__ import annotations
import warnings
from collections.abc import Iterable, Mapping, Sequence
from dataclasses import dataclass, field, replace
from types import MappingProxyType
from typing import Any, Literal, NamedTuple
import numpy as np
from .._inference_status import (
PARAMETER_STATUSES,
InferenceStatus,
is_labelled,
precedent_status,
status_record,
supplies_inference,
)
from .._support import UNCLASSIFIED, Classification, Support, support_for
from .._typing import BoolArray, FloatArray, IntArray
from ..exceptions import refuse_inference
from ..fluctuation.iterative import InitialFit
from ..fluctuation.submodel import Submodel
from ..msm import link_for, solve_projection
from ..utils.bounds import OutcomeScaler, bound
from .cluster import (
cluster_reference_df,
influence_variance,
stacked_second_moment_variance,
)
from .delta import (
Transform,
log_odds_ratio_influence,
log_ratio_influence,
wald_ci,
wald_statistic,
)
__all__ = [
"ClusterReference",
"CorrectionParts",
"CovarianceRule",
"ICParts",
"InferenceStatus",
"ParameterEstimate",
"Scale",
"WaldTest",
"atc_estimate",
"att_estimate",
"bootstrap_column",
"bootstrap_is_diagnostic",
"bootstrap_line",
"counterfactual_mean_parts",
"counterfactual_means",
"has_t_reference",
"ipsi_means",
"is_labelled",
"make_estimate",
"median_estimates",
"minimum_reference_df",
"missing_outcome_correction_parts",
"msm_coefficients",
"natural_course_mean",
"policy_means",
"ratio_estimates",
"reduced_correction_parts",
"reduced_corrections",
"reference_label",
"regime_means",
"restamp_parameter_statuses",
"spread_name",
"stamp_inference",
"supplies_inference",
"support_column",
"support_notes",
]
Scale = Literal["level", "difference", "ratio", "fraction"]
#: How an estimate's joint covariance is read off its influence curve.
CovarianceRule = Literal["centered", "second_moment"]
# ``InferenceStatus`` and ``supplies_inference`` live in the leaf module
# :mod:`cleverly._inference_status`, beside the table of what each status says, and are
# re-exported here, where the documentation names them. ``inference`` is a field with a
# class-level default on :class:`ParameterEstimate`, so a constructor call that omits it
# declares ``"influence_curve"``.
#: The name each spread column takes when the package supplies no inference. The one
#: table of "which number under which name": every frame, record and printed label that
#: publishes a spread of the reported curve reads its name through :func:`spread_name`,
#: and an inferential name with no entry here cannot be published for a diagnostic
#: estimate at all. ``p_value`` has no entry, because a diagnostic estimate reports no
#: p-value under any name.
_DIAGNOSTIC_NAMES: Mapping[str, str] = MappingProxyType(
{
# ParameterEstimate.to_dict, truncation_curve and missingness_tilt.
"std_err": "plugin_std_err",
"ci_lower": "plugin_interval_lower",
"ci_upper": "plugin_interval_upper",
# A percentile interval is a confidence interval too. The numbers stay, under
# the "percentile range" wording ``TMLEResult.summary`` prints for them.
"bootstrap_ci_lower": "bootstrap_range_lower",
"bootstrap_ci_upper": "bootstrap_range_upper",
# The bootstrap spread is the sample standard deviation of the replicate
# estimates, which is a standard error only where the status supplies inference.
# ``bootstrap_std_err`` is the ``to_dict`` key and ``bootstrap se`` the summary label.
"bootstrap_std_err": "bootstrap_sd",
"bootstrap se": "bootstrap sd",
# RefutationTest.to_frame, which names the accessor it read.
"std_error": "plugin_std_error",
# StudyResult.to_frame.
"mean_std_error": "mean_plugin_std_error",
# NuisanceDiagnostics.repeat_spread_frame and its printed table.
"reported_standard_error": "plugin_standard_error",
"ratio_to_standard_error": "ratio_to_plugin_standard_error",
"reported se": "plugin se",
"sd/se": "sd/plugin se",
# The noun the coverage-study verdict traces a shortfall to.
"reported standard error": "plug-in standard error",
# The unit a refutation detail and a score-check verdict count in.
"standard errors": "plug-in standard errors",
"influence-curve standard errors": "plug-in standard errors",
# CVTargeting.to_frame and its printed table.
"cv_std_err": "cv_plugin_std_err",
"pooled_std_err": "pooled_plugin_std_err",
"cv std_err": "cv plugin_std_err",
# SensitivityBounds.to_dict and summary, robustness_value, and the assessment's
# robustness row. The one-sided limits reuse ``ci_lower`` and ``ci_upper``.
"robustness_value_ci": "robustness_value_plugin_interval",
"rva": "rv_plugin_interval",
"confidence-limit value": "plug-in interval value",
"one-sided CIs": "one-sided plug-in intervals",
"confidence bound": "plug-in interval limit",
}
)
def spread_name(name: str, status: InferenceStatus) -> str:
"""The name a spread of the reported curve is published under at an inference status.
An estimate the package supplies inference for keeps every inferential name. A
diagnostic estimate publishes the same number under the name
:data:`_DIAGNOSTIC_NAMES` gives it, which makes no coverage claim.
Parameters
----------
name : str
The inferential name, as an ordinary fit publishes it.
status : str
The inference status of the estimate the number belongs to. One of
:data:`InferenceStatus`, which the
:doc:`inference reference </technical-reference/inference>` lists.
Returns
-------
str
``name`` under ``"influence_curve"``, and its diagnostic name otherwise.
Raises
------
KeyError
When ``status`` is a diagnostic status and ``name`` has no diagnostic name. A
new inferential column has to be entered in the table before a diagnostic fit
can publish it.
"""
if supplies_inference(status):
return name
return _DIAGNOSTIC_NAMES[name]
def bootstrap_column(name: str, status: InferenceStatus, diagnostic: bool) -> str:
"""The name a bootstrap column is published under.
A bootstrap column keeps its inferential name only when the estimate's status supplies
inference *and* the bootstrap is not a diagnostic (:func:`bootstrap_is_diagnostic`).
Otherwise it takes the diagnostic name of :data:`_DIAGNOSTIC_NAMES`.
:func:`spread_name` keeps its one key, the status.
Parameters
----------
name : str
The inferential name, such as ``"bootstrap_ci_lower"`` or ``"bootstrap se"``.
status : str
The inference status of the estimate. One of :data:`InferenceStatus`.
diagnostic : bool
Whether the bootstrap is published as a diagnostic.
Returns
-------
str
``name``, or its diagnostic name.
"""
if supplies_inference(status) and not diagnostic:
return name
return _DIAGNOSTIC_NAMES[name]
#: Why a bootstrap percentile interval is printed as a range on a fit whose status
#: supplies inference: the fit declares the bootstrap a diagnostic.
DIAGNOSTIC_BOOTSTRAP_NOTE = (
"a diagnostic: the collaborative fit selects its propensity model on each resample, and "
"no registered study measures the coverage of that interval"
)
#: The pointwise bootstrap interval: ``"percentile"`` (the ``np.quantile`` limits) or
#: ``"bca"`` (bias-corrected and accelerated, read as order statistics).
BootstrapInterval = Literal["percentile", "bca"]
#: The printed name of each bootstrap interval.
BOOTSTRAP_INTERVAL_LABELS: Mapping[str, str] = MappingProxyType(
{"percentile": "percentile", "bca": "BCa"}
)
def bootstrap_is_diagnostic(estimate: ParameterEstimate) -> bool:
"""Whether the bootstrap of an estimate is published as a diagnostic.
The one rule. A bootstrap is a diagnostic when the estimate's status withholds its
interval, when the fit declares the bootstrap a diagnostic
(:attr:`BootstrapSummary.diagnostic`), or when the fit reads fewer clusters than the
floor (:attr:`BootstrapSummary.below_floor`). The cluster count decides the last case, not
the status, so a fit with estimated weights and few clusters is a diagnostic too. A
bootstrap that none of these makes a diagnostic publishes the inferential names, and its
:attr:`BootstrapSummary.support` says whether a registered row measures it.
Parameters
----------
estimate : ParameterEstimate
An estimate that carries a bootstrap summary.
Returns
-------
bool
``True`` when the summary prints a percentile range and ``to_dict`` publishes the
diagnostic names.
"""
summary = estimate.bootstrap
return not estimate.supplies_inference or (
summary is not None and (summary.diagnostic or summary.below_floor)
)
def bootstrap_line(estimate: ParameterEstimate) -> str:
"""The printed bootstrap spread and interval of one estimate.
One body for both result summaries. The labels come from
:func:`~cleverly.inference.influence.bootstrap_column` under
:func:`bootstrap_is_diagnostic`, so a fit whose status supplies no inference, an
estimate with a t reference, or a bootstrap not licensed as inference prints
``bootstrap sd`` and a range.
Parameters
----------
estimate : ParameterEstimate
An estimate that carries a bootstrap summary.
Returns
-------
str
The line, without the estimate's name.
"""
summary = estimate.bootstrap
assert summary is not None
low, high = summary.ci
label = bootstrap_column("bootstrap se", estimate.inference, bootstrap_is_diagnostic(estimate))
spread = f"{label} {summary.std_error:.4g} "
interval = BOOTSTRAP_INTERVAL_LABELS[summary.interval]
if label == "bootstrap se":
marker = " [experimental]" if summary.support.label == "experimental" else ""
return f"{spread}{interval} CI [{low:.5g}, {high:.5g}]{marker}"
# A bootstrap interval is a confidence interval, and this estimate reports none. The
# same two numbers are printed as a range under the diagnostic framing, rather than
# the bootstrap being refused outright.
if not estimate.supplies_inference:
note = status_record(estimate.inference).bootstrap_note
elif summary.below_floor:
note = status_record("few_cluster_plugin").bootstrap_note
else:
note = DIAGNOSTIC_BOOTSTRAP_NOTE
return f"{spread}{interval} range [{low:.5g}, {high:.5g}] ({note})"
@dataclass(frozen=True)
class WaldTest:
"""A two-sided Wald test of one estimate against a declared null value.
The statistic is ``(inference_value - null_inference_value) / std_error``, with both
values on the estimate's inference scale. That is the log scale for a ratio and the
scale of the transform for a transformed contrast. This is the statistic of
``wald_test(null = )`` in R ``drtmle`` 1.1.2 (``R/test.R``, lines 85 to 183).
Parameters
----------
name : str
The name of the tested estimate.
null : float
The null value on the reported scale, as the caller wrote it.
null_inference_value : float
The null value on the inference scale.
statistic : float
The Wald statistic.
pvalue : float
The two-sided p-value, on the normal reference or on the Student t reference of
the estimate.
scale : {"level", "difference", "ratio", "fraction"}
The reported scale of the estimate.
reference_df : int or None, default=None
The degrees of freedom of the t reference, as
:attr:`ParameterEstimate.reference_df` gives them. ``None`` for the normal
reference.
"""
name: str
null: float
null_inference_value: float
statistic: float
pvalue: float
scale: Scale
reference_df: int | None = None
@dataclass(frozen=True)
class ClusterReference:
r"""The rows each estimate of a clustered report reads, for its Student t reference.
A builder records which rows an estimate reads, and the stamp turns that record into
degrees of freedom. An estimate named in :attr:`stratum_of` reads the rows of its
baseline stratum. Every other estimate reads every row.
An estimate named in :attr:`fold_evaluated` carries the fold-evaluated variance of
:func:`~cleverly.inference.cross_validated_variance`. That variance centres the cluster
totals inside each of the :attr:`validation_folds` folds. Its reference policy uses
:math:`J - V`, the pooled within-group count. Below
:data:`~cleverly._inference_status.FEW_CLUSTER_THRESHOLD` clusters such an estimate takes
:math:`\min(J - 2, J - V)`. The registered studies measure this policy. The count is
not an exact degrees-of-freedom result for heterogeneous cluster variances.
Parameters
----------
cluster : ndarray of int
The cluster label of each row.
weights : ndarray of float or None, default=None
The observation weight of each row, or ``None`` for an unweighted fit.
strata : ndarray of int or None, default=None
The baseline stratum code of each row, or ``None`` for a fit without strata.
stratum_of : mapping of str to int, default={}
The stratum code of each stratum estimate, keyed by estimate name.
fold_evaluated : frozenset of str, default=frozenset()
The names whose variance is the fold-evaluated one.
validation_folds : int, default=0
The number of validation folds that variance averages over.
"""
cluster: IntArray
weights: FloatArray | None = None
strata: IntArray | None = None
stratum_of: Mapping[str, int] = field(default_factory=lambda: MappingProxyType({}))
fold_evaluated: frozenset[str] = frozenset()
validation_folds: int = 0
def with_fold_evaluated(self, names: Iterable[str], validation_folds: int) -> ClusterReference:
r"""This record with ``names`` read as fold-evaluated estimates over that many folds.
Parameters
----------
names : iterable of str
The names whose variance is the fold-evaluated one.
validation_folds : int
The number of validation folds.
Returns
-------
ClusterReference
A copy that applies :math:`\min(J - 2, J - V)` to ``names``.
"""
return replace(
self, fold_evaluated=frozenset(names), validation_folds=int(validation_folds)
)
def reference_df(self, name: str) -> int | None:
"""The degrees of freedom of the estimate called ``name``.
Parameters
----------
name : str
The estimate name.
Returns
-------
int or None
:func:`~cleverly.inference.cluster.cluster_reference_df` over the rows the
estimate reads.
"""
code = self.stratum_of.get(name)
if code is None:
df = cluster_reference_df(self.cluster, self.weights)
else:
assert self.strata is not None
df = cluster_reference_df(self.cluster, self.weights, rows=self.strata == code)
if df is None or name not in self.fold_evaluated:
return df
# ``df`` is J - 2, so J - V is ``df + 2 - V``. The one-cluster-fold refusal keeps
# V at most J / 2, so the value stays positive.
return min(df, df + 2 - self.validation_folds)
def stamp_inference(
estimates: Mapping[str, ParameterEstimate],
status: InferenceStatus,
reference: ClusterReference | None = None,
classification: Classification | None = None,
) -> dict[str, ParameterEstimate]:
"""Give each estimate of a report the status, reference and support its fit gives.
The one stamp. Under a status that withholds its interval, every estimate declares
``status`` and keeps the normal reference of its diagnostic. Under a status that shows its
interval, every estimate that carries no status of its own declares ``status`` and takes
the degrees of freedom of the rows it reads. Under ``"influence_curve"`` every estimate
keeps the status it was built with, and, when ``reference`` is given, takes the degrees of
freedom too. ``TMLE._retarget_detailed`` applies it to the fit's estimates, and
:meth:`~cleverly.estimators.CVTargeting.stamped` to both fold-level reports.
When ``classification`` is given, each estimate also takes the
:class:`~cleverly._support.Support` that :func:`~cleverly._support.support_for` computes
from its final status and reference.
Parameters
----------
estimates : mapping of str to ParameterEstimate
The report, keyed by estimand.
status : str
One of :data:`InferenceStatus`.
reference : ClusterReference or None, default=None
The rows each estimate reads, on a clustered fit. ``None`` keeps every
:attr:`ParameterEstimate.reference_df` as built.
classification : Classification or None, default=None
The classification of the fit. ``None`` keeps every :attr:`ParameterEstimate.support`
as built.
Returns
-------
dict of str to ParameterEstimate
``estimates`` in the same order. Every array and number is unchanged.
"""
stamped: dict[str, ParameterEstimate]
if not supplies_inference(status):
stamped = {
name: replace(estimate, inference=status, reference_df=None)
for name, estimate in estimates.items()
}
else:
stamped = {}
for name, estimate in estimates.items():
changes: dict[str, Any] = {}
if is_labelled(status) and estimate.inference == "influence_curve":
changes["inference"] = status
if reference is not None:
changes["reference_df"] = reference.reference_df(name)
stamped[name] = replace(estimate, **changes) if changes else estimate
if classification is None:
return stamped
return {
name: replace(
estimate,
support=support_for(classification, estimate.inference, estimate.reference_df),
)
for name, estimate in stamped.items()
}
def has_t_reference(estimates: Mapping[str, ParameterEstimate]) -> bool:
"""Whether any estimate of a report carries a Student t reference.
Parameters
----------
estimates : mapping of str to ParameterEstimate
The report, keyed by estimand.
Returns
-------
bool
``True`` when one :attr:`ParameterEstimate.reference_df` is set. A summary then
prints a ``df`` column, and a fit skips its default simultaneous band.
"""
return any(estimate.reference_df is not None for estimate in estimates.values())
def reference_label(estimate: ParameterEstimate) -> str:
"""The ``df`` cell a summary prints for one estimate.
Parameters
----------
estimate : ParameterEstimate
A reported estimate.
Returns
-------
str
The degrees of freedom, or ``"normal"`` for the normal reference.
"""
return "normal" if estimate.reference_df is None else str(estimate.reference_df)
def minimum_reference_df(estimates: Sequence[ParameterEstimate]) -> int | None:
"""The degrees of freedom of an estimate derived from ``estimates``.
The reference policy takes the smallest of the inputs' degrees of freedom, with the
normal reference counting as infinite. When the inputs share a reference, the contrast
retains it. This policy does not establish conservative coverage for correlated inputs
with different references.
Parameters
----------
estimates : sequence of ParameterEstimate
The inputs of the derived estimate.
Returns
-------
int or None
The smallest :attr:`ParameterEstimate.reference_df`, or ``None`` when every input
uses the normal reference.
"""
values = [estimate.reference_df for estimate in estimates if estimate.reference_df is not None]
return min(values) if values else None
[docs]
@dataclass(frozen=True)
class ParameterEstimate:
"""A point estimate with everything needed to do inference on it.
Parameters
----------
name : str
Stable alias of the estimand.
psi : float
The estimate, on the outcome's original scale.
influence_curve : ndarray
Per-observation influence curve on the *inference* scale -- the original
outcome scale for levels and differences, the log scale for ratios.
variance : float
Variance on the inference scale.
n : int
Number of observations behind the estimate.
n_clusters : int
Independent clusters behind it, which is what the variance divides by.
scale : {"level", "difference", "ratio", "fraction"}
Scale the estimate is reported on, which decides how the interval is built.
alpha : float
Significance level of :attr:`ci`.
log_psi : float or None
Present for ratios only: the estimate on the log scale.
bootstrap : BootstrapSummary or None
Bootstrap summary, when ``n_bootstrap > 0`` was requested.
covariance_rule : {"centered", "second_moment"}
How joint covariance and delta-method contrasts are computed from the influence
curve. ``"centered"`` uses the sample covariance at the observation or cluster
unit, which can differ from :attr:`variance` on a fold-evaluated or repeated fit.
``"second_moment"`` uses the raw second moment that the stacked cross-fitted
natural-course mean reports as its variance.
inference : str
Whether the package supplies inference for this estimate. One of
:data:`InferenceStatus`. The
:doc:`inference reference </technical-reference/inference>` lists each status
and its reason.
``"influence_curve"`` is the ordinary case, and :attr:`std_error`, :attr:`ci`
and :attr:`pvalue` answer. A status whose record shows its interval answers the
same three accessors, and the estimate carries the experimental label. At any
other status those three accessors raise with the status's reason, and
:attr:`plugin_std_error` and :attr:`plugin_interval` report the retained
diagnostic.
transform : Transform or None
A monotone map that the interval and the test are computed on, set by
``contrast(..., transform=)``. ``None`` for every estimate a fit reports.
reference_df : int or None
Degrees of freedom of the Student t reference of :attr:`ci`, :attr:`pvalue` and
:attr:`plugin_interval`, or ``None`` for the normal reference. A clustered fit
with fewer than 40 clusters with positive weight mass sets it to the cluster
count minus 2, per estimate.
support : Support
What the package states about the evidence behind this estimate: the label
``"supported"`` or ``"experimental"``, the reason, the evidence and the support row.
A fit sets it from its classification. An estimate built outside a classified fit
is experimental. ``repr`` and ``summary()`` print the word ``experimental`` for an
experimental estimate.
Attributes
----------
std_error : float
plugin_std_error : float
inference_value : float
ci : tuple of float
plugin_interval : tuple of float
pvalue : float
supplies_inference : bool
See Also
--------
cleverly.estimators.TMLEResult : The mapping these are read out of.
cleverly.inference.influence_variance : The variance an estimate is given.
cleverly.inference.simultaneous_bands : Joint bands over several of these.
cleverly.inference.WaldTest : What :meth:`wald_test` returns at a declared null.
Examples
--------
Read an estimate from a fitted result:
>>> from sklearn.linear_model import LinearRegression, LogisticRegression
>>> from cleverly import ATE, CausalStudy, PointTreatment
>>> from cleverly.datasets import make_linear_ate
>>> frame, _ = make_linear_ate(n=80, seed=1)
>>> study = CausalStudy(
... frame,
... design=PointTreatment(
... outcome="Y", treatment="A", adjustment=("W1", "W2", "W3", "W4")
... ),
... )
>>> result = study.identify(ATE()).estimate(
... outcome_learner=LinearRegression(),
... treatment_learner=LogisticRegression(max_iter=1000),
... cross_fit=False,
... random_state=0,
... )
>>> estimate = result["ate"]
>>> estimate.name, estimate.n
('ate', 80)
>>> estimate.ci[0] < estimate.ci[1]
True
Test the estimate against a declared null value:
>>> test = estimate.wald_test(null=1.0)
>>> test.null, test.pvalue < 1.0
(1.0, True)
"""
name: str
psi: float
influence_curve: FloatArray
variance: float
n: int
n_clusters: int
scale: Scale = "difference"
alpha: float = 0.05
log_psi: float | None = None
bootstrap: BootstrapSummary | None = None
covariance_rule: CovarianceRule = "centered"
inference: InferenceStatus = "influence_curve"
transform: Transform | None = None
reference_df: int | None = None
support: Support = UNCLASSIFIED
def _plugin_std_error(self) -> float:
"""The plug-in standard error, under any inference status.
One body for :attr:`std_error` and :attr:`plugin_std_error`, so the retained
diagnostic is the same arithmetic on the same curve rather than a second
implementation of it.
"""
if not np.isfinite(self.variance) or self.variance < 0:
return float("nan")
return float(np.sqrt(self.variance))
def _wald_interval(self) -> tuple[float, float]:
"""The Wald interval of the plug-in curve, under any inference status.
One body for :attr:`ci` and :attr:`plugin_interval`, for the same reason
:meth:`_plugin_std_error` has one.
"""
center = self.inference_value
df = self.reference_df
if self.transform is not None:
low, high = wald_ci(center, self._plugin_std_error(), self.alpha, df=df)
if not (np.isfinite(low) and np.isfinite(high)):
return (float("nan"), float("nan"))
# Sorted, so a decreasing transform still gives an ordered interval.
mapped = sorted(
(float(self.transform.inverse(low)), float(self.transform.inverse(high)))
)
return (mapped[0], mapped[1])
if self.scale == "ratio":
low, high = wald_ci(center, self._plugin_std_error(), self.alpha, df=df)
return (float(np.exp(low)), float(np.exp(high)))
return wald_ci(center, self._plugin_std_error(), self.alpha, df=df)
def _null_inference_value(self, null: float | None) -> tuple[float, float]:
"""The null on the reported scale and on the inference scale."""
with np.errstate(divide="ignore", invalid="ignore"):
return self._map_null(null)
def _map_null(self, null: float | None) -> tuple[float, float]:
if null is None:
null = 1.0 if self.scale == "ratio" else 0.0
if self.transform is not None:
mapped = float(self.transform.forward(null))
if not np.isfinite(mapped):
raise ValueError(
f"the default null {null:g} has no value under the transform "
f"{self.transform.name!r}. Call wald_test(null=...) with a null "
"inside the transform's domain"
)
return float(null), mapped
return float(null), 0.0
null = float(null)
if self.transform is not None:
mapped = float(self.transform.forward(null))
if not np.isfinite(mapped):
raise ValueError(
f"the null {null:g} has no value under the transform "
f"{self.transform.name!r}; pass a null inside the transform's domain"
)
return null, mapped
if self.scale == "ratio":
if not null > 0.0:
raise ValueError(f"a ratio-scale null must be positive; got {null:g}")
return null, float(np.log(null))
return null, null
def _plugin_wald_test(self, null: float | None) -> WaldTest:
"""The Wald test of the plug-in curve, under any inference status.
One body for :meth:`wald_test`, :attr:`pvalue` and the rejection rate
:class:`~cleverly.validation.CoverageStudy` measures, so a study that measures the
diagnostic evaluates the refused property's arithmetic rather than a second rule.
"""
reported, mapped = self._null_inference_value(null)
statistic, pvalue = wald_statistic(
self.inference_value, mapped, self._plugin_std_error(), df=self.reference_df
)
return WaldTest(
name=self.name,
null=reported,
null_inference_value=mapped,
statistic=statistic,
pvalue=pvalue,
scale=self.scale,
reference_df=self.reference_df,
)
def _plugin_pvalue(self) -> float:
"""The two-sided p-value at the scale's default null, under any inference status."""
return self._plugin_wald_test(None).pvalue
def _reported_pvalue(self) -> float:
"""The default test for a report, with a missing test represented by ``nan``."""
try:
return self._plugin_pvalue()
except ValueError:
# A transformed estimate can have a valid interval without a default test.
# Its default null lies outside the transform's domain; wald_test(null=...)
# still lets the caller declare an interior null.
return float("nan")
@property
def supplies_inference(self) -> bool:
"""Whether the package supplies inference for this estimate.
``True`` when :attr:`inference` is ``"influence_curve"`` or a status that shows its
interval, which is when :attr:`std_error`, :attr:`ci` and :attr:`pvalue` answer.
"""
return supplies_inference(self.inference)
[docs]
def spread_columns(self, *, pvalue: bool = True) -> dict[str, float]:
"""The standard error, interval and p-value, each under its status-appropriate name.
An estimate the package supplies inference for returns ``std_err``,
``ci_lower``, ``ci_upper`` and ``p_value``. A diagnostic estimate returns the
same standard error and interval under the names :func:`spread_name` gives them,
and no p-value. The numbers are the bodies :attr:`plugin_std_error` and
:attr:`plugin_interval` read, so they are identical under either status.
Parameters
----------
pvalue : bool, default=True
Whether to include ``p_value`` for an estimate that supplies inference. A
frame with no p-value column passes ``False``.
Returns
-------
dict of str to float
The columns, in the order ``to_dict`` publishes them.
"""
low, high = self._wald_interval()
columns = {
spread_name("std_err", self.inference): self._plugin_std_error(),
spread_name("ci_lower", self.inference): low,
spread_name("ci_upper", self.inference): high,
}
if pvalue and self.supplies_inference:
columns["p_value"] = self._reported_pvalue()
return columns
@property
def std_error(self) -> float:
"""Standard error on the inference scale.
Raises
------
CapabilityError
When :attr:`supplies_inference` is ``False``. Read
:attr:`plugin_std_error` for the retained diagnostic.
"""
refuse_inference(self.inference, operation=".std_error")
return self._plugin_std_error()
@property
def plugin_std_error(self) -> float:
"""Plug-in standard error of the reported curve, at every inference status.
A diagnostic and not a standard error for an estimate whose
:attr:`supplies_inference` is ``False``. On an ordinary estimate it is
:attr:`std_error` under a name that makes no coverage claim.
"""
return self._plugin_std_error()
@property
def inference_value(self) -> float:
"""Return the estimate on the scale used by its influence curve."""
if self.transform is not None:
return float(self.transform.forward(self.psi))
if self.scale == "ratio":
if self.log_psi is None:
raise ValueError("a ratio-scale estimate requires log_psi for inference")
return float(self.log_psi)
return float(self.psi)
@property
def ci(self) -> tuple[float, float]:
"""Confidence interval at level ``1 - alpha``.
Raises
------
CapabilityError
When :attr:`supplies_inference` is ``False``. Read
:attr:`plugin_interval` for the retained diagnostic.
"""
refuse_inference(self.inference, operation=".ci")
return self._wald_interval()
@property
def plugin_interval(self) -> tuple[float, float]:
"""Wald interval of the plug-in curve, at level ``1 - alpha``.
A diagnostic and not a confidence interval. It is the interval :attr:`ci`
builds, read through a name that makes no coverage claim, so a fit whose
inference the package refuses can still report the spread of the curve it
computed. At a non-inferential status the coverage of this interval is
unestablished, for the reason the
:doc:`inference reference </technical-reference/inference>` gives. The
registered selector studies measure it and report it as a diagnostic.
"""
return self._wald_interval()
@property
def pvalue(self) -> float:
"""Two-sided p-value against the null of no effect.
The null is zero for a level, difference or fraction and one for a ratio
(i.e. zero on the log scale). A transformed estimate maps that null through its
transform. The value is ``wald_test().pvalue``, from the same body.
Raises
------
CapabilityError
When :attr:`supplies_inference` is ``False``. This fit reports no
p-value.
"""
refuse_inference(self.inference, operation=".pvalue")
return self._plugin_pvalue()
[docs]
def wald_test(self, null: float | None = None) -> WaldTest:
"""Test the estimate against a declared null value.
The statistic is computed on the inference scale. A ratio's null is mapped to the
log scale, and a transformed estimate's null is mapped through its transform.
Parameters
----------
null : float or None, default=None
The null value on the reported scale. ``None`` takes the scale's default:
zero for a level, difference or fraction, and one for a ratio.
Returns
-------
WaldTest
The statistic and its two-sided p-value. At ``null=None`` its p-value is
:attr:`pvalue`.
Raises
------
CapabilityError
When :attr:`supplies_inference` is ``False``. The fit reports no test.
ValueError
When a ratio-scale null is not positive, or a transformed estimate's null has
no value under its transform.
Examples
--------
>>> import numpy as np
>>> from cleverly.inference import make_estimate
>>> estimate = make_estimate("ey1", 0.4, np.array([0.1, -0.1, 0.2, -0.2]), n=4)
>>> test = estimate.wald_test(null=0.4)
>>> test.statistic, test.pvalue
(0.0, 1.0)
"""
refuse_inference(self.inference, operation=".wald_test()")
return self._plugin_wald_test(null)
@property
def score(self) -> float:
"""Mean of the influence curve.
Targeting is supposed to drive this to zero; :mod:`cleverly.validation.score`
compares it against the standard error to decide whether it did.
"""
return float(np.mean(self.influence_curve))
[docs]
def with_bootstrap(self, summary: BootstrapSummary) -> ParameterEstimate:
"""Return a copy with bootstrap inference attached.
Parameters
----------
summary : BootstrapSummary
Percentile interval and bootstrap standard error to attach.
Returns
-------
ParameterEstimate
A copy carrying that summary. The analytic interval is unchanged.
"""
return replace(self, bootstrap=summary)
[docs]
def to_dict(self) -> dict[str, Any]:
"""Return a JSON-compatible representation.
``support`` is always present, and is ``"supported"`` or ``"experimental"``. An
experimental estimate adds ``support_reason`` and ``support_evidence``, the evidence
joined by ``;``, and an estimate whose support states limits adds ``support_limits``.
``inference`` is present for every status except ``"influence_curve"``, and names it.
An estimate the package supplies inference for, either at ``"influence_curve"`` or at a
status that shows its interval, emits ``std_err``, ``ci_lower``, ``ci_upper`` and
``p_value``. One whose :attr:`supplies_inference` is ``False`` emits
``plugin_std_err``, ``plugin_interval_lower`` and ``plugin_interval_upper`` in their
place, and no ``p_value``. Its bootstrap columns are ``bootstrap_sd``,
``bootstrap_range_lower`` and ``bootstrap_range_upper`` rather than
``bootstrap_std_err``, ``bootstrap_ci_lower`` and ``bootstrap_ci_upper``. The bootstrap
columns take the same diagnostic names when :func:`bootstrap_is_diagnostic` is
``True``. A fit with a bootstrap adds ``bootstrap_support``, the label of the bootstrap
summary. A transformed estimate adds ``transform`` and ``transform_psi``, the estimate
on the transform's scale. An estimate with a Student t reference adds ``reference_df``,
and its bootstrap columns take the diagnostic names. ``std_err`` is present exactly when
:attr:`supplies_inference` is ``True``, so a programmatic consumer has a total test.
Returns
-------
dict
A JSON-compatible mapping of every reported field.
"""
row: dict[str, Any] = {"estimand": self.name, "psi": self.psi}
row["support"] = self.support.label
if self.support.label == "experimental":
row["support_reason"] = self.support.reason
row["support_evidence"] = ";".join(self.support.evidence)
if self.support.limits:
row["support_limits"] = ";".join(self.support.limits)
if is_labelled(self.inference):
row["inference"] = self.inference
row.update(self.spread_columns())
if self.reference_df is not None:
row["reference_df"] = self.reference_df
row["scale"] = self.scale
if self.log_psi is not None:
row["log_psi"] = self.log_psi
if self.transform is not None:
row["transform"] = self.transform.name
row["transform_psi"] = self.inference_value
if self.bootstrap is not None:
diagnostic = bootstrap_is_diagnostic(self)
for key, value in (
("bootstrap_std_err", self.bootstrap.std_error),
("bootstrap_ci_lower", self.bootstrap.ci[0]),
("bootstrap_ci_upper", self.bootstrap.ci[1]),
):
row[bootstrap_column(key, self.inference, diagnostic)] = value
row["bootstrap_support"] = self.bootstrap.support.label
row["bootstrap_interval"] = self.bootstrap.interval
if self.bootstrap.interval == "bca":
row["bootstrap_z0"] = self.bootstrap.z0
row["bootstrap_acceleration"] = self.bootstrap.acceleration
return row
def __repr__(self) -> str: # pragma: no cover - debugging aid
marker = " [experimental]" if self.support.label == "experimental" else ""
if not self.supplies_inference:
low, high = self.plugin_interval
return (
f"{self.name}: {self.psi:.5g} ({status_record(self.inference).summary_label} "
f"{self.plugin_std_error:.4g}, plug-in interval [{low:.5g}, {high:.5g}]; "
"a diagnostic, and no confidence interval or p-value is reported)" + marker
)
low, high = self.ci
reference = "" if self.reference_df is None else f"t({self.reference_df}) "
return (
f"{self.name}: {self.psi:.5g} (se {self.std_error:.4g}, "
f"{reference}{(1 - self.alpha) * 100:g}% CI [{low:.5g}, {high:.5g}], "
f"p={self._reported_pvalue():.3g})" + marker
)
def support_column(estimates: Mapping[str, ParameterEstimate]) -> list[str] | None:
"""The ``support`` cell of each estimate, or ``None`` when every estimate is supported.
A summary adds the column only when an estimate is experimental, so a fully supported
report prints the table it always printed.
Parameters
----------
estimates : mapping of str to ParameterEstimate
The report, keyed by estimand.
Returns
-------
list of str or None
``"supported"`` or ``"experimental"`` per estimate, in report order, or ``None``.
"""
labels: list[str] = [estimate.support.label for estimate in estimates.values()]
return labels if "experimental" in labels else None
def support_notes(estimates: Mapping[str, ParameterEstimate]) -> list[str]:
"""One ``experimental:`` note per distinct reason among the experimental estimates.
Parameters
----------
estimates : mapping of str to ParameterEstimate
The report, keyed by estimand.
Returns
-------
list of str
The notes in first-seen order. A reason that covers only some estimates names them.
The list is empty when every estimate is supported.
"""
reasons: dict[str, list[str]] = {}
for name, estimate in estimates.items():
if estimate.support.label == "experimental":
reasons.setdefault(estimate.support.reason, []).append(name)
notes = []
for reason, names in reasons.items():
covers_all = len(names) == len(estimates)
head = "experimental" if covers_all else f"experimental ({', '.join(names)})"
notes.append(f"{head}: {reason}")
return notes
@dataclass(frozen=True)
class BootstrapSummary:
"""The replicate distribution of one estimand under the full-refit bootstrap.
The fields describe the construction. The inference status of the estimate decides
the names that :meth:`ParameterEstimate.to_dict` and the result summary publish them
under: at a status that supplies no inference, ``std_error`` is published as
``bootstrap_sd`` and ``ci`` as a range.
Parameters
----------
std_error : float
Sample standard deviation of the finite replicate estimates, ``ddof=1``. ``nan``
with fewer than two finite estimates, as every limit below.
ci : tuple of float
Two-sided limits of the interval named by ``interval``.
ci_one_sided_lower : float
Lower one-sided limit.
ci_one_sided_upper : float
Upper one-sided limit.
n_replicates : int
Replicates with a finite estimate.
n_failed : int
Failed full refits. A derived parameter also counts requested draws with missing,
nonfinite, or out-of-domain derived values.
draws : ndarray
The finite replicate estimates.
diagnostic : bool, default=False
Whether the fit declares the bootstrap interval a diagnostic. ``True`` publishes
the columns under their diagnostic names, through :func:`bootstrap_column`, at every
status. Only the collaborative fit sets it.
below_floor : bool, default=False
Whether the fit reads fewer clusters with positive weight mass than the interval floor
(10 for a point fit, 20 for a longitudinal one). ``True`` publishes the columns under
their diagnostic names at every status, because no registered study measures a
bootstrap interval below the floor.
support : Support, default=UNCLASSIFIED
Whether a registered row measures the interval of the fit that produced it. A row
that lists the interval's kind (``bootstrap`` or ``bootstrap_bca``) makes it
``supported``. Every other fit that publishes the interval makes it ``experimental``.
interval : {"percentile", "bca"}, default="percentile"
Which interval ``ci`` and the one-sided limits are: the ``np.quantile`` percentile
limits, or the bias-corrected and accelerated limits of
:func:`~cleverly.inference.bootstrap.bca_interval`.
z0 : float or None, default=None
The BCa bias correction, on a BCa summary.
acceleration : float or None, default=None
The BCa acceleration, from the estimate's own influence curve (summed within clusters
when clusters were resampled), on a BCa summary.
"""
std_error: float
ci: tuple[float, float]
ci_one_sided_lower: float
ci_one_sided_upper: float
n_replicates: int
n_failed: int
draws: FloatArray
diagnostic: bool = False
below_floor: bool = False
support: Support = UNCLASSIFIED
interval: BootstrapInterval = "percentile"
z0: float | None = None
acceleration: float | None = None
@property
def variance(self) -> float:
return float(self.std_error**2)
#: Why the raw second-moment rule refuses clusters. One text for the estimate constructor
#: and the covariance helpers in :mod:`cleverly.inference.results`.
_SECOND_MOMENT_CLUSTER_REFUSAL = (
"the raw second-moment covariance rule is defined for independent rows only, "
"and this result declares clusters"
)
def make_estimate(
name: str,
psi: float,
influence_curve: FloatArray,
*,
n: int,
cluster: IntArray | None = None,
scale: Scale = "difference",
alpha: float = 0.05,
log_psi: float | None = None,
covariance_rule: CovarianceRule = "centered",
inference: InferenceStatus = "influence_curve",
reference_df: int | None = None,
support: Support = UNCLASSIFIED,
) -> ParameterEstimate:
"""Assemble a :class:`ParameterEstimate`, computing its variance under its rule.
This function owns the map from a covariance rule to a stored variance, so an
estimate's :attr:`~ParameterEstimate.variance` and its declared
:attr:`~ParameterEstimate.covariance_rule` cannot disagree.
Parameters
----------
name : str
Stable alias of the estimand.
psi : float
The estimate, on the outcome's original scale.
influence_curve : ndarray
Per-observation influence curve on the inference scale.
n : int
Number of observations behind the estimate.
cluster : ndarray or None, default=None
``(n,)`` cluster codes, or ``None`` for independent rows.
scale : {"level", "difference", "ratio", "fraction"}, default="difference"
Scale the estimate is reported on.
alpha : float, default=0.05
Significance level of the interval.
log_psi : float or None, default=None
The estimate on the log scale, for a ratio only.
covariance_rule : {"centered", "second_moment"}, default="centered"
``"centered"`` stores :func:`~cleverly.inference.influence_variance`.
``"second_moment"`` stores
:func:`~cleverly.inference.cluster.stacked_second_moment_variance` and refuses a
cluster.
inference : str, default="influence_curve"
Whether the package supplies inference for the estimate. One of
:data:`InferenceStatus`, which the
:doc:`inference reference </technical-reference/inference>` lists. The default
is the ordinary case, so a caller that does not pass it builds an inferential
estimate.
reference_df : int or None, default=None
Degrees of freedom of a Student t reference, which a derived estimate computes
with :func:`minimum_reference_df`. ``None`` takes the normal reference. A fitted
report's own estimates take theirs from :func:`stamp_inference`.
support : Support, default=UNCLASSIFIED
The support the estimate carries. A caller that builds an estimate inside a classified
fit passes the support of that fit.
Returns
-------
ParameterEstimate
The estimate, declaring ``covariance_rule``, ``inference``, ``reference_df`` and
``support``.
"""
ic = np.asarray(influence_curve, dtype=float).reshape(-1)
if covariance_rule == "second_moment":
if cluster is not None:
raise ValueError(_SECOND_MOMENT_CLUSTER_REFUSAL)
variance = stacked_second_moment_variance(ic)
elif covariance_rule == "centered":
variance = influence_variance(ic, cluster)
else:
raise ValueError(
f"covariance_rule must be 'centered' or 'second_moment'; got {covariance_rule!r}"
)
n_clusters = n if cluster is None else int(np.unique(cluster).size)
return ParameterEstimate(
name=name,
psi=float(psi),
influence_curve=ic,
variance=variance,
n=n,
n_clusters=n_clusters,
scale=scale,
alpha=alpha,
log_psi=log_psi,
covariance_rule=covariance_rule,
inference=inference,
reference_df=reference_df if supplies_inference(inference) else None,
support=support,
)
def restamp_parameter_statuses(
per_repeat: Sequence[Mapping[str, ParameterEstimate]],
) -> list[dict[str, ParameterEstimate]]:
"""Give every draw of an estimand the precedent of its draws' parameter statuses.
A parameter status (:data:`~cleverly._inference_status.PARAMETER_STATUSES`) belongs to
one estimate of one draw, such as an arm whose stacked curve failed on that draw's
split. The median over draws reads every draw, so it takes the status of the most
restrictive draw. An estimate holds a parameter status only beside an
``"influence_curve"`` fit status (:func:`stamp_inference`), so the restamp applies when
every other draw of the estimand also reports ``"influence_curve"``. Draws with
different fit statuses keep them, and :func:`median_estimates` refuses the
combination: that guard finds a real defect.
The order of the stamps on a repeated fit is this restamp, then the median. A label
that reads the final status, such as a support row, is computed after both.
Parameters
----------
per_repeat : sequence of mappings
One estimate mapping per draw.
Returns
-------
list of dict
The draws, with each restamped estimate replaced. The input is not changed.
"""
out = [dict(report) for report in per_repeat]
if len(out) < 2:
return out
shared = [name for name in out[0] if all(name in report for report in out)]
for name in shared:
statuses = [report[name].inference for report in out]
if not any(status in PARAMETER_STATUSES for status in statuses):
continue
fit_level = {status for status in statuses if status not in PARAMETER_STATUSES}
if not fit_level <= {"influence_curve"}:
continue
status = precedent_status(statuses)
for report in out:
report[name] = replace(report[name], inference=status)
return out
def median_estimates(
per_repeat: Sequence[Mapping[str, ParameterEstimate]],
) -> dict[str, ParameterEstimate]:
r"""Combine repeated cross-fitting estimates by their median.
Each entry is one complete draw of the cross-fitting split. The point estimate is
the median over draws. Ratios use the log scale, which is also their inference scale.
The variance follows Chernozhukov et al. (2018, equation 3.14) and zEpid's repeated
cross-fit TMLE implementation:
.. math::
\widetilde\sigma^2 = \operatorname{median}_r\left\{
\widehat\sigma_r^2 + (\widehat\psi_r - \widetilde\psi)^2\right\}.
This rule includes within-draw uncertainty and between-draw displacement. For an odd
repeat count, the reported marginal influence curve comes from the draw at the median
point. For an even count, it is the mean of the two central draws' curves. The curve
supports marginal diagnostics. The adjusted variance, not the curve's empirical
variance, supplies the interval.
A draw whose estimate holds a parameter status passes it to the median
(:func:`restamp_parameter_statuses`). Draws that declare different fit statuses for one
estimate raise.
Parameters
----------
per_repeat : sequence of mappings
One estimate mapping per draw. A one-element sequence returns its input unchanged.
Returns
-------
dict
Median-combined estimates for names present in every draw.
"""
if not per_repeat:
raise ValueError("median_estimates needs at least one repeat")
if len(per_repeat) == 1:
return dict(per_repeat[0])
per_repeat = restamp_parameter_statuses(per_repeat)
shared = [name for name in per_repeat[0] if all(name in report for report in per_repeat)]
dropped = [name for name in per_repeat[0] if name not in shared]
if dropped:
# Silence here would report a median over a subset of the draws under the same
# name as a median over all of them.
warnings.warn(
f"{', '.join(dropped)} was estimated in some cross-fitting repeats but not "
"in all of them -- a draw whose targeting step failed, or a fold with no "
"units in the conditioning arm -- so it is omitted from the median report "
"rather than combined over fewer draws than the rest.",
UserWarning,
stacklevel=2,
)
out: dict[str, ParameterEstimate] = {}
for name in shared:
parts = [report[name] for report in per_repeat]
rules = {part.covariance_rule for part in parts}
if len(rules) != 1:
raise ValueError(
f"{name} declares different covariance rules {sorted(rules)} across "
"repeats; a median over draws needs one rule"
)
statuses = {part.inference for part in parts}
if len(statuses) != 1:
raise ValueError(
f"{name} declares different inference statuses {sorted(statuses)} across "
"repeats; a median over draws is inferential output or it is not, and "
"combining the two would report one draw's refusal under another's name"
)
supports = {part.support for part in parts}
if len(supports) != 1:
raise ValueError(
f"{name} carries different support records across repeats; a median over "
"draws estimates one configuration, and combining two would label the median "
"with one draw's evidence"
)
scale = parts[0].scale
if scale == "ratio":
# ``points`` and ``median_point`` are on the reporting scale, which for a
# ratio is the log scale. The displacement term below reads them, so a ratio
# gets its between-draw spread in log units and not in ratio units.
points = np.asarray([part.log_psi for part in parts], dtype=float)
median_point = float(np.median(points))
log_psi: float | None = median_point
psi = float(np.exp(median_point))
else:
points = np.asarray([part.psi for part in parts], dtype=float)
median_point = float(np.median(points))
psi = median_point
log_psi = None
order = np.argsort(points, kind="stable")
middle = len(parts) // 2
central = order[middle : middle + 1] if len(parts) % 2 else order[middle - 1 : middle + 1]
influence_curve = np.mean(
[np.asarray(parts[int(index)].influence_curve, dtype=float) for index in central],
axis=0,
)
variance = float(
np.median(
[
part.variance + (point - median_point) ** 2
for part, point in zip(parts, points, strict=True)
]
)
)
out[name] = ParameterEstimate(
name=name,
psi=psi,
influence_curve=np.asarray(influence_curve, dtype=float),
variance=variance,
n=parts[0].n,
n_clusters=parts[0].n_clusters,
scale=scale,
alpha=parts[0].alpha,
log_psi=log_psi,
# The stored variance is the median-adjusted rule above under either
# declaration. The declared rule still decides how covariance reads the curve.
covariance_rule=parts[0].covariance_rule,
inference=parts[0].inference,
reference_df=minimum_reference_df(parts),
support=parts[0].support,
)
return out
def _expect(submodel: Submodel, group: str) -> None:
"""Refuse a submodel from the wrong group.
Every curve here reads a fluctuation built for one score equation, and reading the
wrong one would not raise on its own -- the shapes line up and the arithmetic would
quietly answer for a different parameter. Hence a guard at the top of each, and
hence one function rather than the eight copies of the same two lines it replaces.
"""
if submodel.group != group:
raise ValueError(f"expected the {group!r} submodel; got {submodel.group!r}")
def _residual(outcome: FloatArray, targeted: InitialFit, observed: BoolArray | None) -> FloatArray:
r"""``Delta * (Y - Q*(A, W))`` on the scaled outcome scale."""
y = np.asarray(outcome, dtype=float).reshape(-1)
residual = y - targeted.observed
if observed is None:
return residual
return np.where(np.asarray(observed, dtype=bool), residual, 0.0)
class ArmMean(NamedTuple):
r"""One counterfactual mean :math:`E[Y(a)]` and its influence curve.
Both on the *scaled* outcome scale; :meth:`~cleverly.targets.TargetContext.finish`
maps back.
Also the carrier for the other things a fluctuation produces one of per parameter: a
working model's coefficients (:func:`msm_coefficients`) and the conditional effects
(:func:`att_estimate`). What is shared is the shape -- an estimate and the curve it
is entitled to -- which is all anything downstream reads.
"""
psi: float
influence_curve: FloatArray
def natural_course_mean(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
weights: FloatArray,
observed: BoolArray,
) -> ArmMean:
r"""The MAR natural-course mean and efficient influence curve.
.. math::
\hat\psi = P_n m^*(A,W), \qquad
D^*(O) = \frac{\Delta}{\pi(A,W)}\{Y-m^*(A,W)\}+m^*(A,W)-\psi.
No treatment mechanism appears: the plug-in integrates over the empirical joint
distribution of ``(A, W)`` rather than substituting a fitted treatment law.
"""
_expect(submodel, "natural_course")
w = np.asarray(weights, dtype=float).reshape(-1)
prediction = np.asarray(targeted.observed, dtype=float)
psi = float(np.average(prediction, weights=w))
residual = _residual(outcome, targeted, observed)
return ArmMean(psi, w * (submodel.observed[:, 0] * residual + prediction - psi))
def counterfactual_means(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
weights: FloatArray,
observed: BoolArray | None = None,
corrections: Mapping[float, FloatArray] | None = None,
) -> dict[float, ArmMean]:
r"""Every counterfactual mean and its influence curve, keyed by arm.
.. math::
\hat\Psi_a = \frac1n \sum_i \bar Q^*(a, W_i),
\qquad
D_a^*(O) = h_a(A, W)\,\{Y - \bar Q^*(A, W)\} + \bar Q^*(a, W) - \Psi_a
A mapping rather than the ``(psi1, IC1, psi0, IC0)`` tuple this used to return: the
four-tuple could not describe a third arm, and unpacking it was the last thing in the
estimand layer that counted the arms.
``submodel`` must be the *unweighted* mean submodel even when the fluctuation was fit
in weighted form: the influence curve is defined by the true clever covariate, not by
the reparameterisation used to fit it.
``corrections`` are :func:`reduced_corrections`' :math:`D^*_Q + D^*_g`, **subtracted**,
for a fit that solved the doubly-robust equations. ``None`` -- every fit that did not --
leaves the expression below untouched character for character, which matters: see the
comment on the sum. The subtraction cannot move :math:`\hat\Psi`, since the targeting
drove all three empirical means to zero; what it moves is the variance, which is the
whole of what that variant buys.
"""
_expect(submodel, "mean")
w = np.asarray(weights, dtype=float).reshape(-1)
residual = _residual(outcome, targeted, observed)
out: dict[float, ArmMean] = {}
for arm in targeted.levels:
prediction = targeted.arms[arm]
psi = float(np.average(prediction, weights=w))
# Summed in this association deliberately. Splitting the bracket to reuse
# ICParts here would be mathematically identical and would change the last bit of
# every influence curve, because floating-point addition is not associative. The
# decomposition is a diagnostic (`counterfactual_mean_parts`); the estimation path
# keeps the arithmetic its regression fixtures were built against.
if corrections is None:
out[arm] = ArmMean(psi, w * (submodel.column_for(arm) * residual + prediction - psi))
else:
out[arm] = ArmMean(
psi,
w
* (
submodel.column_for(arm) * residual
+ prediction
- psi
- np.asarray(corrections[arm], dtype=float)
),
)
return out
def reduced_corrections(
outcome: FloatArray,
targeted: InitialFit,
treatment: FloatArray,
reduced: Any,
propensity: FloatArray,
*,
bounds: tuple[float, float],
guard: tuple[str, ...],
) -> dict[float, FloatArray]:
r""":math:`D^*_Q + D^*_g` per arm, the terms doubly-robust inference subtracts.
.. math::
D^*_g &= \frac{Q_r(a, W)}{g^*(a|W)}\,\{1_a - g^*(a|W)\} \\
D^*_Q &= 1_a\,\frac{g_{r,2}(a|W)}{g_{r,1}(a|W)}\,\{Y - \bar Q^*(a, W)\}
The second display is the default univariate form. The bivariate construction uses
:math:`D^*_Q=1_a(g_r-g^*)/(g^*g_r)\{Y-\bar Q^*\}` instead; the mechanism correction is
unchanged.
Randomized missing-outcome fits do not enter this two-correction function. They use
:func:`missing_outcome_correction_parts`, which keeps the paper's treatment,
observation and outcome blocks separate.
In either setting the reported curve is :math:`D^* - D^*_Q - D^*_g`. **Minus**, both
of them: that is
what ``drtmle`` computes and what Theorem 1 derives (see below), and a sum is the
plausible transcription error. Since the
targeting drove all three empirical means to zero, the combination cannot move the point
estimate however the signs go -- it moves only the variance, so nothing that reports
:math:`\hat\Psi` can catch getting this wrong. ``tests/unit/test_influence_drtmle.py``
carries the sign as a negative control at deliberately wrong nuisances, which is the
only place it is visible: at the truth :math:`Q_r` and the outcome-correction numerator
vanish row by row and the reported curve *equals* :math:`D^*` array for array.
**Validity is not efficiency.** Under misspecification the canonical gradient at
:math:`P_0` is still :math:`D^*`; what this subtraction leaves is the **estimator's**
asymptotic influence function at the nuisance limits, and the estimator is generally not
efficient there. So a doubly-robust interval is not a narrower one and not an efficient
one. With one nuisance wrong it is valid only in the cases
``docs/technical-reference/dr-tmle/theorem.md`` lists. When both nuisances are
consistent, :math:`Q_r` and the outcome-correction numerator converge to zero and the
curve approaches the efficient curve. At the true nuisance functions, both corrections
vanish row by row.
Worth saying explicitly because the numbers invite the opposite reading -- in the guide's
own worked example the corrected standard error is the smaller of the two.
**The sign was challenged and the challenge is answered.** The 2016 working paper's §3.1
*display* defines the mechanism-side correction with a leading minus,
:math:`D_A = -Q_r/g\,(A - g)`, and Theorem 1 subtracts :math:`D_A` -- which read together
would make the theorem's curve :math:`D^* + Q_r/g\,(A - g)`, the opposite of this
function's. The sign discrepancy is **resolved in favour of this implementation** on
the paper's own appendices: each derives its block as
:math:`P_0[\text{term}] = -(P_n - P_0)D + P_n D + (\text{second order})`, an identity
satisfiable only with :math:`D` equal to the **positive** term, so appendix A forces
:math:`D_A = +Q_r/g\,(A - g)` and appendix B forces :math:`D_Y` to be the positive
quantity below. Theorem 1's own :math:`\sigma^2_n = P_n\{D^* - D_A - D_Y\}^2` is then
exactly what this function's caller computes. The paper also prints :math:`D_Y` twice
with two signs, which is the other reason not to settle it from a display.
``tests/unit/test_theorem_drtmle.py`` checks the appendix step that fixes the
orientation, and pins these arrays against the theorem's terms;
``docs/technical-reference/dr-tmle/theorem.md`` carries the argument.
The published 2017 article confirms both corrected TMLE constructions; the sign adjudication
still rests on the explicit appendix identity plus exact-law arithmetic rather than on a
second transcription of the same formula.
**And the decomposition is now pinned against a perturbation of the law**, which is the
check every other estimand in this package gets and the one this variant could not have
until a fixture was made wrong on purpose. In each half of the union model -- one
nuisance consistent, the other a declared constant -- with the reduced regressions
saturated, the corrected curve *is* the efficient influence function row for row:
.. math::
\frac{1}{g_1} - \frac{g_{r,2}}{g_{r,1}} = \frac{1}{g_0},
\qquad
D^* (\bar Q^*, g_0) - D^*_g = D^*(\bar Q_0, g_0),
the second because the :math:`\bar Q^*` in the residual and the one inside
:math:`Q_r = \bar Q_0 - \bar Q^*` cancel. The identity holds because the consistent
nuisance separates every cell of ``tests/discrete_law.py``. On a continuous law the
conditioning term of ``docs/technical-reference/dr-tmle/theorem.md`` can remain
(``tests/unit/test_conditioning_remainder_drtmle.py``).
``tests/unit/test_influence_gateaux_drtmle.py`` checks the identity against
``tests/discrete_law.py``'s complex-step derivative -- from a real fit as well as from
longhand, at ``rtol=0``. A flipped sign misses by ``0.55`` to ``2.8`` against a
``1e-12`` window, which was watched rather than argued.
And note exactly what the reference to ``drtmle`` is: **no number this package produces
has been compared against that package's, and none will be.** There is no cross-language
test here or in CI and none is planned. So "what that package computes" is a statement
about a formula read out of its source -- **provenance** -- and not about agreement with
anything it returns. That is a decision and not a gap: both implementations descend from
one source, so agreement would be evidence about the transcription and blind to exactly
the error above, which is why the sign was settled by the appendices instead.
``docs/architecture-invariants.md``, *Validation and evidence*, carries the reasoning for
using independent validation and the condition that would reopen it.
**The two terms are built by** :func:`reduced_correction_parts` **and added here**, so
that :func:`~cleverly.validation.drtmle.correction_check` takes each one's empirical
mean from this expression rather than from a second copy of it. A second copy is how
an identity check comes to agree with a curve neither of them is. A fit reaches the
sum through :meth:`CorrectionParts.total` rather than through this function, which is
where the two terms are also reported apart; this is the sum written down once.
**One correction per equation the fit actually solved**, which is what ``guard`` selects
and is the crossing ``guard=`` has everywhere else: :math:`D^*_g` is equation (9)'s, the
one the ``"Q"`` guard adds, and :math:`D^*_Q` is equation (10)'s, the one ``"g"`` adds.
A fit guarding one nuisance subtracts one term. Subtracting both would leave the unsolved
equation's arbitrary mean in the curve, measured at
:math:`2.8\times10^{-3}` on a ``guard=("g",)`` fit against a :math:`7.7\times10^{-6}` bar
with **no** row clipped, so it is not the bounded-mechanism centring defect in another
guise. The derivation was
already in the repository -- ``tests/unit/test_remainder_drtmle.py`` adds each correction
only under the guard whose equation removes it, and shows that two guards over-correct on
an exact law; ``tests/unit/test_drtmle_fit.py`` fits one end to end. ``DRTMLE``'s own
default is both guards, so the ordinary fit is unaffected.
Parameters
----------
targeted, propensity:
The **targeted** regression and the targeted :math:`g^*(a_1 \mid W)` as an ``(n,)``
array, both as the alternation left them. Reading the initial mechanism here would
report a curve for a fit nobody ran, since the equations were solved at the tilted
one.
reduced:
The :class:`~cleverly.estimators.reduced.ReducedSet` the equations were finally
solved against -- the refit, not the fit's own. Typed loosely to keep this module
free of :mod:`cleverly.estimators`, which imports it.
bounds:
The same mechanism truncation the clever covariates divided by.
guard:
Which of the two extra equations this fit solved, in ``DRTMLE.guard``'s vocabulary,
and so which corrections belong in its curve. **Required, with no default**, which
is the point: an earlier caller did not pass this, and a default of both would
make that caller's mistake the fallback for the next one.
"""
return reduced_correction_parts(
outcome,
targeted,
treatment,
reduced,
propensity,
bounds=bounds,
guard=guard,
).total()
@dataclass(frozen=True)
class CorrectionParts:
r"""The two corrections kept apart, plus what the mechanism truncation absorbed.
:func:`reduced_corrections` is the sum of the first two and is what the reported curve
subtracts. They are built here rather than there so that
:func:`~cleverly.validation.drtmle.correction_check` can take each one's empirical mean
**from the same expression the curve carries** -- an identity checked against a second
implementation of the same formula is not an identity, and this is the one class of
defect that check exists to catch.
Attributes
----------
d_g, d_q:
Rowwise :math:`D^*_g(a)` and :math:`D^*_Q(a)` per arm, on the ``[0, 1]`` scaled
outcome that :math:`Q_r` and the fluctuation's residual both live on.
clip_bias:
Rowwise :math:`Q_r(a, W)/g^b(a|W)\,\{g(a|W) - g^b(a|W)\}` per arm -- the quantity
called :math:`B_{clip}`, in the orientation that expression defines it in.
It is **exactly** the difference between the mechanism score the
alternation solves, at the raw tilted :math:`g^*`, and the mean of the
:math:`D^*_g` above, which truncates :math:`g^*` in its residual as well as in its
denominator -- *negated*, since the residual is
:math:`1_a - g` in one and :math:`1_a - g^b` in the other:
.. math::
P_n[w\,D^*_g] - S_g^{\text{stored}} = P_n[w\,B_{clip}]
Zero on every row the bound leaves alone, so its mean is zero whenever nothing
clips. It is a diagnostic and not a correction: nothing subtracts it.
**It is now zero on every fit**, and that is the intended convention rather than a broken
column. :func:`~cleverly.fluctuation.mechanism.solve_bounded_mechanism` runs at the
``DRTMLE``
call sites, so the mechanism the alternation carries forward is already truncated
and there is no second array left for this to measure the distance to. What still
varies between draws, and what a fixture is chosen on, is
:attr:`~cleverly.validation.drtmle.CorrectionRow.margin`.
clipped:
Which rows the mechanism truncation binds on. Empty on every converged bounded fit for the
same reason: a converged bounded tilt lies inside the bounds. On record because
"the identity holds" was uninformative on a draw where the bound never bit -- the
degeneracy that hid the centring mismatch -- and because zero here is now the
assertion that the bounded convention is active.
guard:
Which equations the fit solved, and so which of the two terms :meth:`total` puts in
the curve. **Both arrays are built whatever it says**, because the term a fit does
*not* subtract is exactly what
:func:`~cleverly.validation.drtmle.correction_check` reports as a diagnostic -- and
it has to come from this expression rather than from a second copy of it, for the
reason above. It travels here rather than being read off the fluctuation twice, so
that :func:`~cleverly.estimators.tmle.correction_parts` is the one place the guard
is copied off the record and the curve and the check cannot select differently.
"""
d_g: dict[float, FloatArray]
d_q: dict[float, FloatArray]
clip_bias: dict[float, FloatArray]
clipped: BoolArray
guard: tuple[str, ...]
d_a: dict[float, FloatArray] | None = None
d_m: dict[float, FloatArray] | None = None
d_y: dict[float, FloatArray] | None = None
def total(self) -> dict[float, FloatArray]:
"""What the curve subtracts: one correction per equation :attr:`guard` solved.
:math:`D^*_Q + D^*_g` under both guards, and under both alone -- the ``"Q"`` guard
is what adds equation (9), so it is what puts :math:`D^*_g` here, and ``"g"`` adds
equation (10) and puts :math:`D^*_Q` here. Membership rather than equality, since
``DRTMLE`` validates the guard's contents and not its order.
The both-guards line below is character for character the single expression this
replaced. Re-associating it -- through ``sum()``, or a comprehension over a list
of terms -- moves the last bit of every ordinary doubly-robust curve, for the
reason :func:`counterfactual_means` records about its own arithmetic.
An empty guard **raises**. Such a fit fits no reduced regressions at all and never
reaches here; returning zeros would make it the plain estimator recovered by a
branch, which is the thing ``DRTMLE._nuisances``' short circuit exists to avoid.
"""
has_q, has_g = "Q" in self.guard, "g" in self.guard
# Aliased rather than copied on the single-guard branches: nothing mutates what
# this returns -- the curve subtracts it and `_slice_fit`'s sibling indexing copies.
if has_q and has_g:
return {arm: np.asarray(self.d_g[arm] + self.d_q[arm], dtype=float) for arm in self.d_g}
if has_q:
return {arm: np.asarray(self.d_g[arm], dtype=float) for arm in self.d_g}
if has_g:
return {arm: np.asarray(self.d_q[arm], dtype=float) for arm in self.d_g}
raise ValueError(
"a doubly-robust curve needs at least one guard; an empty guard fits no "
"reduced regressions and must not reach the corrections at all"
)
def reduced_correction_parts(
outcome: FloatArray,
targeted: InitialFit,
treatment: FloatArray,
reduced: Any,
propensity: FloatArray,
*,
bounds: tuple[float, float],
guard: tuple[str, ...],
) -> CorrectionParts:
""":func:`reduced_corrections`' two terms before they are added, and the clipping bias.
Every argument means what it means there. Both terms are built whatever ``guard``
says; it is :meth:`~cleverly.inference.influence.CorrectionParts.total` that selects,
in the same association the single expression used, so the reported curve is unchanged
to the last bit.
**Complete-data fits only**, and there is no missing-outcome mask to pass: a fit with
``delta=`` and a non-empty guard builds the five reductions and goes to
:func:`missing_outcome_correction_parts` instead, while an empty guard fits no
reductions and forms no corrections at all.
"""
y = np.asarray(outcome, dtype=float).reshape(-1)
a = np.asarray(treatment, dtype=float).reshape(-1)
raw = np.asarray(propensity, dtype=float)
# The carrier rule: one column is the two-arm complement form. A composite mechanism
# arrives as ``(n, K)`` at every arm count, because its columns do not sum to one.
if raw.ndim == 1:
if len(reduced.arms) != 2:
raise ValueError(f"a one-column mechanism requires two arms; got {list(reduced.arms)}")
raw1 = raw.reshape(-1)
g1 = bound(raw1, float(bounds[0]), float(bounds[1]))
mechanism = {reduced.arms[0]: 1.0 - g1, reduced.arms[1]: g1}
# The complement rather than a separately clipped array, exactly as
# `Propensity.bounded` and `reduced_mechanism_covariate` take it.
untruncated = {reduced.arms[0]: 1.0 - raw1, reduced.arms[1]: raw1}
clipped = np.asarray(raw1 != g1, dtype=bool)
else:
if raw.shape != (y.size, len(reduced.arms)):
raise ValueError(
f"the targeted mechanism must be ({y.size}, {len(reduced.arms)}); got {raw.shape}"
)
bounded = bound(raw, float(bounds[0]), float(bounds[1]))
mechanism = {arm: bounded[:, j] for j, arm in enumerate(reduced.arms)}
untruncated = {arm: raw[:, j] for j, arm in enumerate(reduced.arms)}
# Reduced **over the arms**, so this stays one bit per row as the binary branch's
# is and as `clipped` is documented and reported: the arms are clipped column by
# column here, and counting the cells would report "3800 row(s) of 2000" on a
# three-armed fit whose bound binds on two arms at most rows.
clipped = np.asarray((raw != bounded).any(axis=1), dtype=bool)
gr = reduced.bounded_gr1(bounds)
d_g: dict[float, FloatArray] = {}
d_q: dict[float, FloatArray] = {}
clip_bias: dict[float, FloatArray] = {}
for j, arm in enumerate(reduced.arms):
# The residual and the mechanism have to denote the same event -- `mechanism` is
# P(A = a | W), so the indicator is 1(A = a) and nothing else multiplies it. A
# mask on only one side is the canonical-source discrepancy that gated this
# feature, and it is the reason `observed=` is not a parameter here.
indicator = (a == float(arm)).astype(float)
qr = np.asarray(reduced.qr, dtype=float)[:, j]
d_g[arm] = qr / mechanism[arm] * (indicator - mechanism[arm])
# The outcome residual is at the arm this row took, so the indicator already puts
# it at `arm` and nothing further selects rows.
if reduced.reduction == "univariate":
ratio = np.asarray(reduced.gr2, dtype=float)[:, j] / gr[:, j]
else:
# van der Laan's bivariate correction and pinned R `eval_Dstar_Q`:
# I(A=a) / gr * (gr-g) / g * (Y-Q). `mechanism` is the same bounded
# targeted margin used by the ordinary and Qr corrections.
ratio = (gr[:, j] - mechanism[arm]) / (mechanism[arm] * gr[:, j])
d_q[arm] = indicator * ratio * (y - targeted.observed)
clip_bias[arm] = qr / mechanism[arm] * (untruncated[arm] - mechanism[arm])
return CorrectionParts(d_g, d_q, clip_bias, clipped, tuple(guard))
def missing_outcome_correction_parts(
outcome: FloatArray,
targeted: InitialFit,
treatment: FloatArray,
observed: BoolArray,
reduced: Any,
propensity: FloatArray,
missingness: FloatArray,
*,
g_bounds: tuple[float, float],
missingness_bound: float,
guard: tuple[str, ...],
) -> CorrectionParts:
r"""The separate treatment, observation and outcome corrections in the paper.
``propensity`` is the targeted treatment mechanism in the shape the fit tilted it. At two
arms it is ``g_1`` of the one shared tilt, and arm 0 takes the complement. Above two arms
it is the ``(n, K)`` armwise mechanism, and each column is bounded separately.
"""
y = np.asarray(outcome, dtype=float).reshape(-1)
a = np.asarray(treatment, dtype=float).reshape(-1)
delta = np.asarray(observed, dtype=float).reshape(-1)
if len(reduced.arms) == 2:
raw_g = np.asarray(propensity, dtype=float).reshape(-1)
bounded_upper = bound(raw_g, float(g_bounds[0]), float(g_bounds[1]))
g_a = np.column_stack([1.0 - bounded_upper, bounded_upper])
clipped_g = np.asarray(raw_g != bounded_upper, dtype=bool)
else:
raw_k = np.asarray(propensity, dtype=float)
if raw_k.shape != (y.size, len(reduced.arms)):
raise ValueError(
f"the targeted mechanism must be ({y.size}, {len(reduced.arms)}); got {raw_k.shape}"
)
g_a = bound(raw_k, float(g_bounds[0]), float(g_bounds[1]))
# Reduced over the arms, so that `clipped` stays one bit per row, as in
# `reduced_correction_parts`.
clipped_g = np.asarray((raw_k != g_a).any(axis=1), dtype=bool)
raw_m = np.asarray(missingness, dtype=float)
g_m = bound(raw_m, float(missingness_bound), 1.0)
gamma_a = reduced.bounded_gamma_a(g_bounds)
gamma_m = reduced.bounded_gamma_m(missingness_bound)
w2 = np.asarray(reduced.r_a, dtype=float) / (gamma_a * gamma_m)
w2 += np.asarray(reduced.r_m, dtype=float) / gamma_m
e = np.asarray(reduced.e, dtype=float)
d_a: dict[float, FloatArray] = {}
d_m: dict[float, FloatArray] = {}
d_y: dict[float, FloatArray] = {}
d_g: dict[float, FloatArray] = {}
zeros: dict[float, FloatArray] = {}
for j, arm in enumerate(reduced.arms):
indicator = (a == float(arm)).astype(float)
d_a[arm] = e[:, j] / g_a[:, j] * (indicator - g_a[:, j])
d_m[arm] = indicator * e[:, j] / (g_a[:, j] * g_m[:, j]) * (delta - g_m[:, j])
d_y[arm] = indicator * delta * w2[:, j] * (y - np.asarray(targeted.observed, dtype=float))
d_g[arm] = np.asarray(d_a[arm] + d_m[arm], dtype=float)
zeros[arm] = np.zeros_like(d_g[arm])
clipped = np.asarray(clipped_g | np.any(raw_m != g_m, axis=1), dtype=bool)
return CorrectionParts(
d_g=d_g,
d_q=d_y,
clip_bias=zeros,
clipped=clipped,
guard=tuple(guard),
d_a=d_a,
d_m=d_m,
d_y=d_y,
)
def policy_means(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
weights: FloatArray,
observed: BoolArray | None = None,
mixing: FloatArray | None = None,
rows: BoolArray | None = None,
terms: FloatArray | None = None,
) -> dict[float, ArmMean]:
r"""Each policy's mean and influence curve, keyed by policy code.
.. math::
\hat\Psi_r = \frac1n \sum_i \bar Q^*\bigl(d_r(A_i, W_i), W_i\bigr),
\qquad
D_r^*(O) = h_r(A, W)\,\{Y - \bar Q^*(A, W)\}
+ \bar Q^*\bigl(d_r(A,W), W\bigr) - \Psi_r
**Why this is not** :func:`regime_means` **at the induced density.** A shift
:math:`d` induces the density :math:`g^d(b \mid w) = \sum_{a : d(a,w)=b} g(a \mid w)`,
and the two parameters are equal:
.. math::
E[\bar Q(d(A,W), W)]
= E_W\Bigl[\sum_a g(a \mid W) \bar Q(d(a,W), W)\Bigr]
= E_W\Bigl[\sum_b g^d(b \mid W) \bar Q(b, W)\Bigr].
The clever covariates are equal too, entry for entry. The **influence curves are
not**. A regime's plug-in term is :math:`\sum_b g^d(b \mid W)\,\bar Q^*(b, W)`, a
function of :math:`W` alone; a shift's is :math:`\bar Q^*(d(A,W), W)`, which reads the
dose the unit actually received. The two agree only in conditional expectation given
:math:`W`, so their difference
.. math::
\bar Q^*(d(A,W), W) - E\bigl[\bar Q^*(d(A,W), W) \mid W\bigr]
is mean zero given :math:`W` and uncorrelated with the regime curve -- the residual
half of that curve is centred given :math:`(A, W)`, and its plug-in half is
:math:`W`-measurable. Hence the exact identity
.. math::
\operatorname{Var}(D^*_{\text{mtp}})
= \operatorname{Var}(D^*_{\text{regime}})
+ \operatorname{Var}\bigl(\bar Q(d(A,W),W) - E[\bar Q(d(A,W),W) \mid W]\bigr).
An MTP has variance at least as large as the known stochastic regime with the same mean.
The added term is a conditional variance, and it is positive wherever the shift moves the
dose -- the price of an intervention that reads the natural value of treatment. So this
must not delegate to :func:`regime_means`, and
``tests/unit/test_influence_gateaux_shift.py`` keeps a negative control that fails if
someone later makes it.
Arithmetically this *is* :func:`counterfactual_means` with the arm replaced by the
shift: ``arm_columns`` is populated for the ``mtp`` group, so ``column_for`` answers,
and ``targeted.arms[r]`` already holds :math:`\bar Q^*` at the shifted dose. The
bracket below is summed in the same association for the same reason the comment there
gives.
**Data-tuned caps (X40 part (e)).** ``rows`` restricts component ``c`` to the rows
``rows[:, c]``: a fold-average policy's component ``k`` is the map fitted on fold
``k``'s training rows, and its mean :math:`\hat\psi_k` is the weighted mean over fold
``k``'s held-out rows. Row ``i`` of fold ``k`` carries
:math:`(W / W_k) w_i (D_i - \hat\psi_k)`, with :math:`W` the total and :math:`W_k`
the fold's weight mass, so the policy's ``mixing`` row of :math:`1/V` gives Hubbard,
Kherad-Pajouh and van der Laan's (2016) Theorem 1 variance. ``terms`` adds, to policy
``s``'s curve, the chain-rule term :math:`\partial\Phi_n' \mathrm{IF}` of a
population target (``(n, S)``, zero for every other policy).
"""
_expect(submodel, "mtp")
w = np.asarray(weights, dtype=float).reshape(-1)
residual = _residual(outcome, targeted, observed)
out: dict[float, ArmMean] = {}
for index, code in enumerate(targeted.levels):
prediction = targeted.arms[code]
if rows is None or bool(np.all(rows[:, index])):
psi = float(np.average(prediction, weights=w))
curve = w * (submodel.column_for(code) * residual + prediction - psi)
else:
held = np.asarray(rows[:, index], dtype=bool)
psi = float(np.average(prediction[held], weights=w[held]))
scale = float(np.sum(w)) / float(np.sum(w[held]))
curve = np.where(
held,
scale * w * (submodel.column_for(code) * residual + prediction - psi),
0.0,
)
out[code] = ArmMean(psi, np.asarray(curve, dtype=float))
if mixing is not None:
out = _mixed_policy_means(out, mixing)
if terms is None:
return out
return {
code: ArmMean(
part.psi,
np.asarray(part.influence_curve, dtype=float) + np.asarray(terms[:, int(code)]),
)
for code, part in out.items()
}
def _mixed_policy_means(out: dict[float, ArmMean], mixing: FloatArray) -> dict[float, ArmMean]:
"""Each policy's mean and curve as the known-weight sum of its components'."""
# A randomized policy's mean is the known-weight sum of its components' means, and so
# is its curve: the randomizer law does not depend on P, so integrating it out is
# linearity of expectation and adds no term (Diaz et al. 2023, Section 2).
components = [out[float(index)] for index in range(len(out))]
mixed: dict[float, ArmMean] = {}
for row, probabilities in enumerate(np.asarray(mixing, dtype=float)):
used = [(float(p), components[c]) for c, p in enumerate(probabilities) if p > 0.0]
psi = sum(p * part.psi for p, part in used)
curve = sum(p * np.asarray(part.influence_curve, dtype=float) for p, part in used)
mixed[float(row)] = ArmMean(float(psi), np.asarray(curve, dtype=float))
return mixed
def regime_means(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
regimes: FloatArray,
weights: FloatArray,
observed: BoolArray | None = None,
) -> dict[float, ArmMean]:
r"""Every regime's counterfactual mean and influence curve, keyed by regime code.
.. math::
\hat\Psi_r = \frac1n \sum_i \sum_a g^\star_r(a \mid W_i)\, \bar Q^*(a, W_i),
\qquad
D_r^*(O) = h_r(A, W)\,\{Y - \bar Q^*(A, W)\}
+ \sum_a g^\star_r(a \mid W)\,\bar Q^*(a, W) - \Psi_r
The counterpart of :func:`counterfactual_means` for the ``regime`` fluctuation, and
the same influence curve with :math:`\mathbb 1\{a = v\}` in place of
:math:`g^\star_r` -- which is why a static regime reproduces the arm result.
Two details are not interchangeable with the arm version. The mixture over arms is
taken **after** fluctuation, on :math:`\bar Q^*` rather than :math:`\bar Q^0`: a
plug-in of the *targeted* distribution is what makes this a substitution estimator
that solves the score equation. And the residual column is read off the submodel
directly, because :attr:`~cleverly.fluctuation.submodel.Submodel.arm_columns` is
empty here -- no column belongs to an arm, so
:meth:`~cleverly.fluctuation.submodel.Submodel.column_for` has nothing to answer.
``regimes`` is the ``(n, K, R)`` density; ``submodel`` must be the *unweighted*
``regime`` submodel even when the fluctuation was fit in weighted form, for the
reason :func:`counterfactual_means` gives.
"""
_expect(submodel, "regime")
star = np.asarray(regimes, dtype=float)
w = np.asarray(weights, dtype=float).reshape(-1)
residual = _residual(outcome, targeted, observed)
levels = targeted.levels
if star.shape[1] != len(levels):
raise ValueError(
f"regimes has {star.shape[1]} arm column(s) but the targeted fit has "
f"{len(levels)} arm(s) {list(levels)}"
)
if star.shape[2] != submodel.dim:
raise ValueError(
f"regimes describes {star.shape[2]} regime(s) but the submodel has "
f"{submodel.dim} column(s)"
)
# (n, K) of the targeted predictions, columns in the same order the densities use.
predictions = np.column_stack([targeted.arms[level] for level in levels])
out: dict[float, ArmMean] = {}
for index in range(star.shape[2]):
mixture = np.einsum("ij,ij->i", star[:, :, index], predictions)
psi = float(np.average(mixture, weights=w))
out[float(index)] = ArmMean(
psi, w * (submodel.observed[:, index] * residual + mixture - psi)
)
return out
def ipsi_means(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
incremental: Any,
treatment: FloatArray,
weights: FloatArray,
observed: BoolArray | None = None,
) -> dict[float, ArmMean]:
r"""Every tilt's counterfactual mean and influence curve, keyed by tilt code.
.. math::
\hat\Psi_r = \frac1n \sum_i m_r(W_i),
\qquad
m_r(W) = \sum_a q_{r}(a \mid W)\, \bar Q^*(a, W)
D_r^*(O) = h_r(A, W)\,\{Y - \bar Q^*(A, W)\}
+ h_r(A, W)\,\{\bar Q^*(A, W) - m_r(W)\}
+ m_r(W) - \Psi_r
with :math:`q_r(a) = \delta_{r,a} g(a)/D_r` and :math:`h_r(a) = \delta_{r,a}/D_r`. The
middle term is the mechanism term: its coefficient is
:math:`\partial m_r/\partial g_b = h_r(b)\{\bar Q^*(b) - m_r\}`, evaluated on the
mechanism score :math:`1\{A=b\} - g_b` and summed over :math:`b`. For two levels it is
:math:`\delta_r\{\bar Q^*(1,W) - \bar Q^*(0,W)\}(A - g^*)/D_r^2`, and that expression is
the one computed there, so the binary arrays are bit-identical to the ones the binary
path produced before the categorical form existed.
**The middle term is the whole reason this is not** :func:`regime_means`. A
:class:`~cleverly.interventions.Stochastic` regime evaluated at the very same density
:math:`q_{r}` has the same mean and, entry for entry, the same clever covariate; its
influence curve is this one without the middle term. Because :math:`q_r` is built out of
:math:`g`, the parameter moves when the mechanism does, and that dependence is a pathwise
derivative the regime curve has no term for.
The gap is not a wash. The extra term is mean zero given :math:`W` and orthogonal to
both halves of the regime curve -- the residual half is centred given :math:`(A, W)`,
the plug-in half is :math:`W`-measurable -- so
.. math::
\operatorname{Var}(D^*_{\text{ipsi}})
= \operatorname{Var}(D^*_{\text{regime}})
+ \operatorname{Var}\!\left(h_r(A)\{\bar Q(A) - m_r\}\right),
an exact decomposition. Treating an incremental intervention as the regime that
induces it therefore does not merely report a different quantity: it reports a
standard error that is too *small*, always.
``tests/unit/test_influence_gateaux_ipsi.py`` keeps that identity as a negative
control, on the same terms ``tests/unit/test_influence_gateaux_shift.py`` does for
the shift axis.
``incremental`` is the :class:`~cleverly.interventions.IPSISet` **as targeted** -- its
:attr:`~cleverly.interventions.IPSISet.propensity` must be the fluctuated mechanism,
not the initial one, or the middle term would be evaluated at a :math:`g` the plug-in
did not use. Carrying the mechanism on the same object as the density is what makes
that impossible to get wrong rather than merely documented.
A tilt of a **known** mechanism (:attr:`~cleverly.interventions.IPSISet.mechanism_known`)
omits the middle term. With :math:`g_0` known the model has no mechanism scores, so the
efficient curve is the full-model curve projected off that tangent space, which removes
exactly the mechanism term. The result is the curve of the known stochastic
regime :math:`q_\delta(g_0)`.
"""
_expect(submodel, "ipsi")
density = np.asarray(incremental.values, dtype=float)
w = np.asarray(weights, dtype=float).reshape(-1)
residual = _residual(outcome, targeted, observed)
levels = targeted.levels
if density.shape[1] != len(levels):
raise ValueError(
f"the tilted density has {density.shape[1]} arm column(s) but the targeted fit "
f"has {len(levels)} arm(s) {list(levels)}"
)
if density.shape[2] != submodel.dim:
raise ValueError(
f"the tilt describes {density.shape[2]} intervention(s) but the submodel has "
f"{submodel.dim} column(s)"
)
# (n, K) of the targeted predictions, columns in the same order the densities use.
predictions = np.column_stack([targeted.arms[level] for level in levels])
a = np.asarray(treatment, dtype=float).reshape(-1)
# Two levels keep the expression the binary path has always used, so that its arrays are
# unchanged to the last bit.
term = _ipsi_term_binary if len(levels) == 2 else _ipsi_term_general
out: dict[float, ArmMean] = {}
for index in range(density.shape[2]):
mixture = np.einsum("ij,ij->i", density[:, :, index], predictions)
psi = float(np.average(mixture, weights=w))
mechanism_term = (
0.0
if incremental.mechanism_known
else term(incremental, predictions, mixture, a, index)
)
curve = submodel.observed[:, index] * residual + mechanism_term + mixture - psi
out[float(index)] = ArmMean(psi, w * curve)
return out
def _ipsi_term_binary(
incremental: Any, predictions: FloatArray, mixture: FloatArray, a: FloatArray, index: int
) -> FloatArray:
r""":math:`\delta\{\bar Q^*(1,W) - \bar Q^*(0,W)\}(A - g^*)/D^2`, Kennedy's mechanism term."""
del mixture
blip = predictions[:, 1] - predictions[:, 0]
treated_residual = a - np.asarray(incremental.propensity, dtype=float)[:, 1]
return np.asarray(incremental.binary_derivative[:, index] * blip * treated_residual)
def _ipsi_term_general(
incremental: Any, predictions: FloatArray, mixture: FloatArray, a: FloatArray, index: int
) -> FloatArray:
r""":math:`h(A)\{\bar Q^*(A, W) - m(W)\}`: the mechanism term for any number of levels."""
at_observed = predictions[np.arange(a.size), np.rint(a).astype(int)]
# The tilt factor at the received level. Not `submodel.observed`, which carries the
# missingness weight of the outcome half and has no place in this term.
return np.asarray(incremental.observed(a)[:, index] * (at_observed - mixture))
def _raw_predictions(
targeted: InitialFit, levels: tuple[float, ...], scaler: OutcomeScaler
) -> FloatArray:
"""``(n, K)`` targeted predictions on the *original* outcome scale, arms in order."""
return scaler.unscale_levels(np.column_stack([targeted.arms[level] for level in levels]))
def msm_coefficients(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
design: FloatArray,
model_weights: FloatArray,
weights: FloatArray,
scaler: OutcomeScaler,
observed: BoolArray | None = None,
link: str = "identity",
) -> dict[float, ArmMean]:
r"""A working model's coefficients and their influence curves, keyed by coefficient.
:math:`\hat\beta` solves the weighted least-squares equation
:math:`U(\beta, \bar Q^*) = 0` of :func:`~cleverly.msm.solve_projection`, which for
the identity link is the closed form
:math:`M^{-1} P_n[\sum_a h\,\varphi\,\bar Q^*]`, and
.. math::
D_\beta^*(O) = M^{-1}\Big[
H(A, W)\,\{Y - \bar Q^*(A, W)\}
+ \sum_a h(a, V)\,\frac{dm}{d\eta}\,\varphi(a, V)\,
\big\{\bar Q^*(a, W) - m(a, V; \hat\beta)\big\}
\Big]
-- the standard M-estimation form, with :math:`M = -\partial U/\partial\beta` the
matrix :func:`~cleverly.msm.solve_projection` returns. :math:`M` carries no *further*
term because :math:`U(\beta_0, \bar Q_0) = 0` at the truth, so the variation of
:math:`M^{-1}` multiplies zero; that argument is untouched by the link, but what
:math:`M` **contains** is not -- under a non-identity link it carries a curvature term
that no saturated working model can exercise. It carries no nuisance either, because
:math:`h` and :math:`\varphi` are known functions; see :mod:`cleverly.msm` for the
estimand and for why the weights have to be known for that to hold.
Both halves are still zero in the sample after targeting, and still for two different
reasons -- but under a link the first is zero only because the fluctuation was solved
at *this* :math:`\hat\beta`, which is what
:func:`~cleverly.estimators.targeting.solve_with_projection` alternates to achieve.
The first term is zero because the ``msm`` fluctuation solved its own score equation.
The second is zero *by construction*: :math:`\hat\beta` is the weighted least-squares
solution against :math:`\bar Q^*`, so the residuals it leaves are orthogonal to
:math:`h\,(dm/d\eta)\,\varphi`. With the identity link the first needs no iteration
either, since the covariate does not mention :math:`\beta` -- which is what makes that
case a one-fluctuation TMLE.
**The projection is solved on the original outcome scale**, unlike every other
estimand here, which works on the scaled outcome and maps back afterwards. A
coefficient vector has no single :class:`Scale` to map back *with*: writing
:math:`\bar Q^*_{\text{raw}} = \ell + r\,\bar Q^*_{\text{scaled}}` gives
:math:`\beta_{\text{raw}} = \ell M^{-1} P_n[\sum_a h \varphi] + r\,
\beta_{\text{scaled}}`, and the first term collapses to "the intercept picks up
:math:`\ell`, the slopes pick up nothing" *only* when the design happens to contain an
intercept column. Solving where the coefficients are reported removes the question
rather than requiring a design to promise something. Nothing is lost: the residual
rescales by the same :math:`r`, so a score that is zero on one scale is zero on the
other. The consequence for callers is that these estimates must **not** go through
:meth:`~cleverly.targets.TargetContext.finish`, which would unscale a second time.
Parameters
----------
design:
``(n, K, p)`` array :math:`\varphi(a, V)`, arms in ``targeted.levels`` order.
model_weights:
``(n, K)`` array :math:`h(a, V)`.
weights:
The ``(n,)`` *observation* weights, which are a different thing: they tilt the
population the projection is taken over, while ``model_weights`` says how the
arms are traded off within it.
scaler:
Maps the targeted predictions back off ``[0, 1]``.
link:
The working model's declared link, by name. A *string* rather than a
:class:`~cleverly.msm.Link`, for the same reason the arrays are passed plainly:
this layer is written against what a saved result carries.
Both arrays are passed plainly rather than as a
:class:`~cleverly.msm.MSMSet`, on the same terms as ``regimes`` and ``policies``: the
inference layer is written against arrays so that it does not depend on the objects
that produced them.
"""
_expect(submodel, "msm")
phi = np.asarray(design, dtype=float)
h = np.asarray(model_weights, dtype=float)
w = np.asarray(weights, dtype=float).reshape(-1)
levels = targeted.levels
if phi.ndim != 3 or phi.shape[1] != len(levels):
raise ValueError(
f"the working model's design must have shape (n, {len(levels)}, p) for arms "
f"{list(levels)}; got {phi.shape}"
)
if phi.shape[2] != submodel.dim:
raise ValueError(
f"the working model has {phi.shape[2]} term(s) but the submodel has "
f"{submodel.dim} column(s)"
)
if h.shape != phi.shape[:2]:
raise ValueError(
f"the working model's weights must have shape {phi.shape[:2]}; got {h.shape}"
)
residual = scaler.unscale_influence(_residual(outcome, targeted, observed))
predictions = _raw_predictions(targeted, levels, scaler)
spec = link_for(str(link))
fit = solve_projection(phi, h, predictions, w, str(link))
beta = fit.beta
fitted = np.asarray(spec.inverse(np.einsum("ijp,p->ij", phi, beta)), dtype=float)
# (n, K, p) -- h * dm/dbeta, which is h * phi under the identity link and so leaves
# that path on exactly the arithmetic it was on before links existed.
weighted_design = phi * (h * np.asarray(spec.slope(fitted), dtype=float))[:, :, None]
# (n, p): the plug-in half, sum_a h * (dm/dbeta) * (Qbar* - m(a, V; beta)).
plugin = np.einsum("ijp,ij->ip", weighted_design, predictions - fitted)
contribution = w[:, None] * (submodel.observed * residual[:, None] + plugin)
influence = np.linalg.solve(fit.jacobian, contribution.T).T
return {
float(j): ArmMean(float(beta[j]), np.ascontiguousarray(influence[:, j]))
for j in range(phi.shape[2])
}
class ICParts(NamedTuple):
r"""The two halves of an influence curve, kept apart.
.. math::
D^*(O) = \underbrace{H(A, W)\,\{Y - Q^*(A, W)\}}_{\text{residual}}
+ \underbrace{Q^*(a, W) - \Psi}_{\text{plug-in}}
They answer different questions, and summing them immediately -- which this code
used to do -- throws that away. A heavy tail in the *residual* term is a
positivity artefact: one unit with a large inverse-propensity weight dominating
the estimating equation. A heavy tail in the *plug-in* term is genuine outcome
heterogeneity, which no amount of truncation will fix. Their relative size also
says how much work targeting is doing.
"""
residual: FloatArray
plugin: FloatArray
#: The doubly-robust correction, ``-(D*_Q + D*_g)``, for a fit that solved the extra
#: score equations -- or the one term of it a single-guard fit solved for, since this
#: is :meth:`CorrectionParts.total` negated and that selects. **Not** the ``guard=``
#: tuple, which shares only the name; zeros for every other fit, which is what keeps
#: ``total`` and
#: ``shares()`` reading exactly as they did before that variant existed. It belongs to
#: neither half above -- it is neither a positivity artefact nor outcome heterogeneity
#: but the term the extra score equations remove from the curve -- and leaving it out
#: would make this decomposition disagree with the curve by the whole of what the
#: variant does.
guard: FloatArray | None = None
@property
def total(self) -> FloatArray:
summed = np.asarray(self.residual + self.plugin, dtype=float)
if self.guard is None:
return summed
return np.asarray(summed + self.guard, dtype=float)
def shares(self) -> dict[str, float]:
"""Each term's share of the influence curve's variance."""
total = float(np.var(self.total))
if total <= 0:
out = {"residual": float("nan"), "plugin": float("nan")}
return out if self.guard is None else {**out, "guard": float("nan")}
shares = {
"residual": float(np.var(self.residual)) / total,
"plugin": float(np.var(self.plugin)) / total,
}
if self.guard is None:
return shares
return {**shares, "guard": float(np.var(self.guard)) / total}
def counterfactual_mean_parts(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
weights: FloatArray,
observed: BoolArray | None = None,
corrections: Mapping[float, FloatArray] | None = None,
) -> dict[float, ICParts]:
"""The decomposed influence curves behind :func:`counterfactual_means`, per arm.
``parts.total`` agrees with the summed curve to floating-point rounding rather
than bit-for-bit: the sum there is bracketed differently, and addition is not
associative. The gap is a few ULP and is asserted in
``tests/unit/test_ic_parts.py``; use :func:`counterfactual_means` for the
estimate and this for the diagnostic.
``corrections`` is what :func:`counterfactual_means` subtracts, and must be passed
whenever it was: this is meant to decompose *the* curve, and a decomposition missing a
term the curve has is worse than no decomposition at all.
"""
_expect(submodel, "mean")
w = np.asarray(weights, dtype=float).reshape(-1)
residual = _residual(outcome, targeted, observed)
return {
arm: ICParts(
w * submodel.column_for(arm) * residual,
w * (targeted.arms[arm] - float(np.average(targeted.arms[arm], weights=w))),
None if corrections is None else -w * np.asarray(corrections[arm], dtype=float),
)
for arm in targeted.levels
}
def att_estimate(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
treatment: FloatArray,
weights: FloatArray,
observed: BoolArray | None = None,
*,
reference: float = 0.0,
) -> dict[float, ArmMean]:
"""The effect among each non-reference arm's own units, on the scaled outcome scale.
One entry per non-reference arm, keyed by that arm -- a single entry, the classic
ATT, when the treatment is binary.
"""
_expect(submodel, "att")
return _conditional_effects(
outcome, targeted, submodel, treatment, weights, observed, reference, conditions_on_arm=True
)
def atc_estimate(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
treatment: FloatArray,
weights: FloatArray,
observed: BoolArray | None = None,
*,
reference: float = 0.0,
) -> dict[float, ArmMean]:
"""The mirror image: each contrast among the *reference* arm's units."""
_expect(submodel, "atc")
return _conditional_effects(
outcome,
targeted,
submodel,
treatment,
weights,
observed,
reference,
conditions_on_arm=False,
)
def _conditional_effects(
outcome: FloatArray,
targeted: InitialFit,
submodel: Submodel,
treatment: FloatArray,
weights: FloatArray,
observed: BoolArray | None,
reference: float,
*,
conditions_on_arm: bool,
) -> dict[float, ArmMean]:
"""``E[Y^a - Y^r | A = c]`` per non-reference arm, with ``c`` the group's population.
``conditions_on_arm`` picks which: the contrast arm itself for the ATT, the reference
arm for the ATC, where every parameter conditions on the same population.
The columns are read by :meth:`~cleverly.fluctuation.submodel.Submodel.contrast_column_for`
rather than positionally, so an implementation that lined the arms up by order would
fail on the three-armed law rather than pass on a coincidence.
"""
a = np.asarray(treatment, dtype=float).reshape(-1)
w = np.asarray(weights, dtype=float).reshape(-1)
residual = _residual(outcome, targeted, observed)
out: dict[float, ArmMean] = {}
for arm in sorted(submodel.contrast_columns):
conditioning = arm if conditions_on_arm else reference
indicator = np.asarray(a == conditioning, dtype=float)
share = float(np.average(indicator, weights=w))
if share <= 0:
raise ValueError(f"no observations in arm {conditioning:g}: the estimand is undefined")
contrast = targeted.arms[arm] - targeted.arms[reference]
psi = float(np.average(contrast, weights=w * indicator))
column = submodel.contrast_column_for(arm)
out[arm] = ArmMean(psi, w * (column * residual + (indicator / share) * (contrast - psi)))
return out
def ratio_estimates(
psi_one: float,
ic_one: FloatArray,
psi_zero: float,
ic_zero: FloatArray,
*,
n: int,
cluster: IntArray | None = None,
alpha: float = 0.05,
which: tuple[str, ...] = ("rr", "or"),
) -> dict[str, ParameterEstimate]:
"""Risk ratio and odds ratio, with inference on the log scale.
Only meaningful for a binary outcome, where the counterfactual means are risks
in ``(0, 1)``; the caller is responsible for not asking otherwise.
"""
out: dict[str, ParameterEstimate] = {}
if "rr" in which:
log_psi, ic = log_ratio_influence(psi_one, ic_one, psi_zero, ic_zero)
out["rr"] = make_estimate(
"rr",
float(np.exp(log_psi)),
ic,
n=n,
cluster=cluster,
scale="ratio",
alpha=alpha,
log_psi=log_psi,
)
if "or" in which:
log_psi, ic = log_odds_ratio_influence(psi_one, ic_one, psi_zero, ic_zero)
out["or"] = make_estimate(
"or",
float(np.exp(log_psi)),
ic,
n=n,
cluster=cluster,
scale="ratio",
alpha=alpha,
log_psi=log_psi,
)
return out
def unscale(
psi: float, ic: FloatArray, scaler: OutcomeScaler, scale: Scale
) -> tuple[float, FloatArray]:
"""Map an estimate and its influence curve back to the original outcome scale.
A level picks up the location shift, a difference does not, and an influence
curve never does -- it is already centred.
"""
if scale == "ratio":
if not scaler.is_identity:
raise ValueError("ratio estimands are only defined for an unscaled (binary) outcome")
return psi, np.asarray(ic, dtype=float)
if scaler.is_identity:
return psi, np.asarray(ic, dtype=float)
value = scaler.unscale_level(psi) if scale == "level" else scaler.unscale_difference(psi)
return value, scaler.unscale_influence(ic)