"""Shared estimator infrastructure: configuration, results and reporting.
Every TMLE variant in this library produces the same shape of result -- a set of
:class:`~cleverly.inference.ParameterEstimate` objects plus the nuisance fits and
targeting details that produced them -- so sensitivity analysis, validation and
reporting are written once, here, rather than per estimator.
"""
from __future__ import annotations
from collections.abc import Callable, Iterator, Mapping, Sequence
from dataclasses import dataclass, field, replace
from functools import cached_property
from typing import TYPE_CHECKING, Any, ClassVar
import numpy as np
from .._inference_status import (
NO_SIMULTANEOUS_BANDS,
NO_T_REFERENCE_BANDS,
T_REFERENCE_NOTE,
InferenceStatus,
is_labelled,
status_record,
supplies_inference,
)
from .._refusal import NO_ALTERNATIVE, Refusal
from .._support import Classification, SupportKey
from .._typing import FloatArray, ParameterAxis
from ..data.causal_data import CausalData
from ..data.treatment_covariates import (
drtmle_restriction_line,
treatment_mechanism_restriction,
)
from ..exceptions import (
CapabilityError,
RefusalReason,
refuse_after_repeats,
refuse_inference,
)
from ..fluctuation.iterative import Fluctuation
from ..inference.bootstrap import BootstrapResult
from ..inference.cluster import (
cluster_sizes,
cluster_weight_mass,
fewest_clusters,
positive_mass_clause,
unequal_cluster_sizes,
)
from ..inference.delta import Transform
from ..inference.influence import (
ClusterReference,
ParameterEstimate,
Scale,
bootstrap_line,
has_t_reference,
reference_label,
spread_name,
stamp_inference,
support_column,
support_notes,
)
from ..inference.multiplier import SimultaneousBands
from ..inference.results import (
RatioKind,
estimate_covariance,
estimate_curves,
ratio_contrast,
ratio_function,
refuse_odds_ratio,
reported_status,
select_estimates,
smooth_contrast,
sole_estimate,
with_derived_bootstrap,
)
from ..learners.crossfit import CrossFitPlan, SplitPlan
from ..protocol import ProtocolDetail, protocol_summary_lines
from ..provenance import Provenance
from ..targets import TARGETS, all_names, resolve_estimands
from ..targets.base import arm_alias, parameter_name
from ..utils.frames import emit_frame
from ..utils.memos import without_memos
from ..utils.text import format_pvalue, format_table
from ._nuisance import NuisanceEstimates, RepeatFit
from .direct_effect import describe as describe_direct_effect
from .targeting import TargetingSpec
if TYPE_CHECKING: # pragma: no cover - typing only
from ..assessment import (
AssessmentReport,
DiagnosticsFacade,
Replayability,
SensitivityFacade,
ValidationReport,
)
from ..validation.score import ScoreCheck
from .ctmle import CTMLEDerivativeFit, CTMLEOutcomeAdaptiveFit, CTMLESelection
__all__ = [
"ALL_ESTIMANDS",
"DEFAULT_ESTIMANDS",
"MEAN_GROUP_ESTIMANDS",
"CVTargeting",
"TMLEConfig",
"TMLEResult",
"TMLEResultSet",
"resolve_estimands",
]
#: Estimands available from a classic point-treatment fit, in report order. Derived
#: from the registry rather than declared here, so a registered target appears in it
#: automatically and the two cannot disagree.
ALL_ESTIMANDS: tuple[str, ...] = all_names()
#: Estimands produced by the two-column ``mean`` fluctuation.
MEAN_GROUP_ESTIMANDS: frozenset[str] = frozenset(
name for name, target in TARGETS.items() if target.group == "mean"
)
#: The default report for a continuous outcome. A binary outcome additionally gets
#: ``rr`` and ``or``, which are only defined when the means are probabilities.
#:
#: Spans every parameter axis, so it is a *listing* rather than a report any one fit
#: produces: no fit reports ``ate`` and ``ey_regime`` and ``ey_policy`` together, because
#: declaring ``interventions=`` or ``policies=`` switches which axis is in scope. Use
#: :func:`~cleverly.targets.default_names` with the fit's axis for the report itself.
DEFAULT_ESTIMANDS: tuple[str, ...] = tuple(
name for name, target in TARGETS.items() if target.in_default_set
)
@dataclass(frozen=True)
class TMLEConfig:
"""A snapshot of the settings a fit actually used.
Recorded on the result so a reported number can always be traced back to the
bounds, folds and submodel that produced it -- including values the estimator
resolved itself, such as ``g_bounds="auto"``.
"""
family: str
#: Everything the targeting step needs that is not data. Held as one object
#: rather than as loose fields so that re-solving a fluctuation -- the truncation
#: curve, the MNAR tilt, the C-TMLE search -- needs the config and not the live
#: estimator that produced it.
targeting_spec: TargetingSpec
targeting_scheme: str
cross_fit: bool
n_folds: int
g_bounds: tuple[float, float]
g_bounds_conditional: tuple[float, float]
missingness_bound: float
q_bounds: tuple[float, float] | None
screen_treatment: bool
estimands: tuple[str, ...]
alpha_sig: float
random_state: int | None = None
n_bootstrap: int = 0
cv_evaluation: bool = False
auto_bounds_n: float | None = None
#: Which mechanisms ``missingness_bound`` was actually applied to -- empty when the
#: fit had neither missing outcomes nor an intermediate variable, in which case the
#: bound exists on the config but never touched anything.
bounded_mechanisms: tuple[str, ...] = ()
#: Whether a treatment mechanism was estimated. ``False`` only for the
#: missing-outcome natural-course mean, whose clever covariate divides by the
#: response probability alone. Recorded rather than re-derived from the estimand
#: list, because "reports ey_obs and bounded some mechanism" is true of fits that
#: did estimate and truncate ``g`` -- and those fits must still print the bound they
#: applied.
fits_treatment: bool = True
#: ``"known"`` when the fit divided by a mechanism the data declares known
#: (:attr:`~cleverly.data.CausalData.known_treatment`), else ``"estimated"``. Recorded
#: so a result read back from disk says which one its interval rests on.
treatment_mechanism: str = "estimated"
#: Whether the outcome regression was fitted once per treatment arm
#: (``outcome_by_arm=True``) rather than once on the covariates and the treatment.
outcome_by_arm: bool = False
#: Whether the treatment is a continuous dose, whose policy ratio reads a density that
#: no bound truncates. A categorical treatment's policy ratio reads the bounded
#: propensity, so its summary prints the bound.
continuous_treatment: bool = False
#: The arm code every contrast estimand is taken against. Part of the *estimand*
#: rather than a setting: ``ate[medium vs low]`` and ``ate[medium vs high]`` are
#: different parameters, so which one was reported has to be recorded alongside the
#: number. ``0.0`` -- the lowest arm -- unless the caller chose otherwise.
reference_arm: float = 0.0
#: What this fit's parameters are indexed by -- see
#: :attr:`~cleverly.Target.parameter_axis`. Recorded rather than inferred from the
#: estimand names, so a result read back from disk can say what it reported without
#: parsing them.
parameter_axis: ParameterAxis = "arm"
#: The fold policy the caller *declared*, as against ``n_folds`` above, which is the
#: count the fit actually ran. The two come apart without leaving a trace otherwise:
#: ``resolve_n_folds`` caps the count at the rarest stratum and again at the cluster
#: count, each with a warning that is gone by the time anyone reads the result.
#: Defaulted so that every existing construction and ``dataclasses.replace`` keeps
#: working untouched.
crossfit: CrossFitPlan = field(default_factory=CrossFitPlan)
#: The support key the fit was classified under, read before any learner ran. Saved
#: with the result, so a result read back from disk says what its label rested on.
#: ``None`` on a config built outside a classified fit.
support_key: SupportKey | None = None
#: The group-A conditions the fit met when it was preflighted. Saved with the result, so a
#: retarget of the fit's nuisances labels its estimates with the conditions of the fit
#: instead of preflighting again on data that may have changed.
support_conditions: tuple[Refusal, ...] = ()
# Read-through to the spec, so the settings appear once and cannot drift.
@property
def fluctuation(self) -> str:
return self.targeting_spec.fluctuation
@property
def targeting(self) -> str:
return self.targeting_spec.targeting
@property
def alpha(self) -> float:
return self.targeting_spec.alpha
@property
def target_weights(self) -> bool:
return self.targeting_spec.target_weights
@property
def estimator_name(self) -> str:
"""Which of the three cross-fitting constructions this fit actually ran.
Worth naming rather than reporting the raw settings: one common update over the
validation losses is the source CV-TMLE targeting step, while a separate update
per fold is an optional extension with different finite-sample behavior.
"""
if not self.cross_fit:
return "TMLE (in-sample nuisances)"
if self.targeting_scheme == "pooled":
return "fold-evaluated CV-TMLE" if self.cv_evaluation else "stacked CV-TMLE (Levy)"
return (
"fold-specific fold-evaluated TMLE"
if self.cv_evaluation
else "fold-specific targeted TMLE"
)
def describe(self) -> list[str]:
"""Human-readable lines for :meth:`TMLEResult.summary`."""
lines = [
f"outcome family: {self.family}; fluctuation: {self.fluctuation} ({self.targeting})",
]
if self.cross_fit:
lines.append(
f"{self.estimator_name}: nuisances cross-fitted over {self.n_folds} folds; "
f"targeting: {self.targeting_scheme}"
)
if self.crossfit.repeated:
# Always shown when it applies, unlike the capped-folds line below: a
# repeated fit reports a different estimator from an ordinary one, and a
# reader who cannot tell them apart cannot compare two summaries.
lines.append(
f" (median over {self.crossfit.repeats} independent draws of the "
"split; variance includes within-draw uncertainty and split dispersion)"
)
else:
lines.append(f"{self.estimator_name}: nuisances fitted in-sample (cross_fit=False)")
if self.cross_fit and self.crossfit.n_folds != self.n_folds:
# Only when they disagree, because agreeing is the ordinary case and a line
# that always appears is a line nobody reads. When they do disagree the
# warning that explained it was emitted at fit time and is long gone, so this
# is the only place a reader can find out that the split was capped.
lines.append(
f" ({self.crossfit.n_folds} folds were declared; the split was capped "
f"at {self.n_folds} by the rarest stratum or the cluster count)"
)
# A policy fit on a dose reads a conditional density, not a propensity: nothing is
# truncated into g_bounds and reporting the bound would name a step that did not
# happen. What bounds a density ratio there is the cap the analyst declared,
# which is part of the estimand and so appears in the parameter names instead.
if not (self.parameter_axis == "policy" and self.continuous_treatment) and (
self.fits_treatment
):
bounds = f"propensity truncated to [{self.g_bounds[0]:.4g}, {self.g_bounds[1]:.4g}]"
if self.auto_bounds_n is not None:
# Named because it is a deliberate divergence from R's rule, and because a
# reader comparing two fits needs to know the bound moved with the weights.
bounds += f" (auto, resolved at the effective n of {self.auto_bounds_n:.0f})"
if self.g_bounds_conditional != self.g_bounds:
bounds += (
f"; ATT/ATC to [{self.g_bounds_conditional[0]:.4g}, "
f"{self.g_bounds_conditional[1]:.4g}]"
)
lines.append(bounds)
if self.treatment_mechanism == "known":
lines.append(
"treatment mechanism: known (declared); no truncation bound moved a known value"
)
if self.bounded_mechanisms:
# The propensity is not the only denominator in the clever covariate, and a
# reader comparing two fits needs to see every bound that shaped the estimate
# -- not just the one with a familiar name.
named = " and ".join(self.bounded_mechanisms)
lines.append(f"{named} truncated to [{self.missingness_bound:.4g}, 1]")
if self.q_bounds is not None:
lines.append(
f"outcome scaled from [{self.q_bounds[0]:.4g}, {self.q_bounds[1]:.4g}] to [0, 1]"
)
if self.target_weights:
lines.append("fluctuation in weighted form (target_weights=True)")
return lines
@dataclass(frozen=True)
class CVTargeting:
"""Fold-level detail from a cross-validated TMLE fit.
The original CV-TMLE evaluates updated fold-specific distributions one fold at a
time. Its usual pooled targeting step fits one common ``epsilon`` by minimising the
equal average of their validation losses; :attr:`epsilon` records that update. Levy's
easy implementation instead stacks the out-of-fold predictions, targets and evaluates
them just like one ordinary TMLE. That stacked construction is the package default,
defined here by Levy (2018). The pinned ``tmle3`` source snapshot in the references
implements the same path by requesting its ``"validation"`` likelihood in
``tmle3_Update``; it is corroborating implementation evidence, not the specification.
The optional ``targeting_scheme="fold"`` extension fits one update per fold; those
coefficients are kept separately in :attr:`fold_epsilon`.
The two reports here are genuinely different estimators, not two views of one.
:attr:`pooled` stitches the updated validation predictions together and evaluates the
estimand once over the whole sample; :attr:`fold_evaluated` evaluates it inside each
fold and averages, which is the construction Zheng & van der Laan (2011) analyse. They
coincide exactly for estimands linear in the targeted predictions at equal fold
sizes -- ``ate``, ``ey1``, ``ey0`` -- and diverge for ``rr``, ``or``, ``att`` and
``atc``, where a ratio of means is not a mean of ratios and the pooled conditional
effects weight by the whole sample's arm share rather than each fold's. Which one
``result[name]`` carries is set by ``TMLE(cv_evaluation=...)``.
The fit builds this object only under ``cv_evaluation=True`` or
``targeting_scheme="fold"``. A default stacked fit returns ``result.cv_targeting``
as ``None``, and ``result[name]`` is then the stacked (Levy) report. :attr:`pooled`
is the whole-sample plug-in at the fitted fluctuation. Under ``cv_evaluation=True``
that fluctuation minimises the fold-reweighted loss, which multiplies each fold's
weights by ``n / (V * sum(w_fold))``. :attr:`pooled` therefore equals the stacked
report only when every fold has equal weight mass, which for unweighted rows means
equal fold sizes. Under ``targeting_scheme="fold"`` it is the stitched report of the
fold-targeted predictions.
Under ``repeats=R`` the fields divide by what they are. The three that are
*estimates* -- :attr:`pooled`, :attr:`fold_evaluated` and :attr:`variance` -- follow every
draw, exactly as the headline report does. The four that are indexed *by fold* --
:attr:`n_folds`, :attr:`fold_sizes`, :attr:`fold_estimates`, :attr:`fold_epsilon` --
describe the first draw alone, because fold 3 of one draw is not fold 3 of another
and there is no correspondence along which to combine them. :attr:`repeats` says how
many draws the first three cover.
Attributes
----------
variance : dict of str to float
The cross-validated variance of Zheng & van der Laan (2011) -- the fold-averaged
second moment of the fold-specific influence curves -- per estimand. Its square
root is the standard error attached to :attr:`fold_evaluated`; the pooled report
carries the ordinary influence-curve one. The two agree when the folds are
balanced and the score equation is solved, so a gap between them is itself
informative. Over ``R`` draws, the median aggregation adds each point's squared
displacement before taking the median across draws. :attr:`std_error` returns
the square root at ``"influence_curve"`` and refuses at every other status;
:attr:`plugin_std_error` returns it at every status.
std_error : dict of str to float
fold_estimates:
Per-estimand tuple of fold-specific plug-in estimates, from the first draw.
Estimands that some fold could not evaluate (no units in the conditioning arm, a
boundary counterfactual mean) are absent, and their omission was warned about at
fit time.
epsilon:
Per-targeting-group common fluctuation coefficient. For the fold-specific
extension this is only the mass-weighted summary reported by ``Fluctuation``.
fold_epsilon:
Per-targeting-group tuple of fold-specific fluctuation coefficients, from the
first draw. Empty for common pooled-validation targeting.
pooled, fold_evaluated:
The two reports, per estimand, combined by their median over every draw.
repeats:
How many cross-fitting draws :attr:`pooled`, :attr:`fold_evaluated` and
:attr:`variance` cover.
inference : str
plugin_std_error : dict of str to float
"""
n_folds: int
fold_sizes: tuple[int, ...]
variance: dict[str, float]
fold_estimates: dict[str, tuple[float, ...]]
epsilon: dict[str, tuple[float, ...]]
fold_epsilon: dict[str, tuple[tuple[float, ...], ...]]
pooled: dict[str, ParameterEstimate] = field(default_factory=dict)
fold_evaluated: dict[str, ParameterEstimate] = field(default_factory=dict)
repeats: int = 1
#: Name of the dataframe backend the fit's data arrived in, as on every other
#: report here.
backend: str | None = None
@property
def inference(self) -> InferenceStatus:
"""The inference status both reports declare.
Read off the stamped estimates of :attr:`pooled` and :attr:`fold_evaluated` rather than
stored, so the fold-level report cannot claim a status its own estimates do not
carry. ``"influence_curve"`` when neither report holds an estimate. A mix raises
:class:`ValueError`, as :attr:`TMLEResult.inference_status` does.
"""
reports = {
f"{label}[{name}]": estimate
for label, report in (
("pooled", self.pooled),
("fold_evaluated", self.fold_evaluated),
)
for name, estimate in report.items()
}
return reported_status(reports)
@property
def std_error(self) -> dict[str, float]:
"""Cross-validated standard error per estimand.
Raises
------
CapabilityError
When :attr:`inference` supplies no inference. Read :attr:`plugin_std_error`
for the same numbers as a diagnostic.
"""
refuse_inference(self.inference, operation="CVTargeting.std_error")
return self.plugin_std_error
@property
def plugin_std_error(self) -> dict[str, float]:
"""The square root of :attr:`variance` per estimand, at every inference status.
:attr:`std_error` under a name that makes no coverage claim. At a
non-inferential status it is a diagnostic and not a standard error.
"""
return {name: float(np.sqrt(value)) for name, value in self.variance.items()}
def stamped(
self,
status: InferenceStatus,
reference: ClusterReference | None = None,
classification: Classification | None = None,
) -> CVTargeting:
"""This report with both fold-level reports declaring ``status``.
The one stamp for the fold-level reports. ``TMLE._retarget_detailed`` applies it
when it builds the report.
Parameters
----------
status : str
One of :data:`~cleverly.inference.influence.InferenceStatus`.
reference : ClusterReference or None, default=None
The rows each estimate reads on a clustered fit, which set the Student t
reference of an inferential estimate. The fold-evaluated report reads it with
:meth:`~cleverly.inference.influence.ClusterReference.with_fold_evaluated`.
classification : Classification or None, default=None
The classification of the fit, which gives each estimate its support.
Returns
-------
CVTargeting
A copy whose :attr:`pooled` and :attr:`fold_evaluated` estimates are stamped by
:func:`~cleverly.inference.influence.stamp_inference`, with every number
unchanged.
"""
fold_reference = (
None
if reference is None
else reference.with_fold_evaluated(self.fold_evaluated, self.n_folds)
)
return replace(
self,
pooled=stamp_inference(self.pooled, status, reference, classification),
fold_evaluated=stamp_inference(
self.fold_evaluated, status, fold_reference, classification
),
)
def to_frame(self, data: CausalData | None = None) -> Any:
"""One row per estimand: both reports, the CV standard error and the spread.
The two standard-error columns take their names from
:func:`~cleverly.inference.influence.spread_name`. At a labelled status the frame
also carries an ``inference`` column naming it. At a status that withholds its
interval the columns read ``cv_plugin_std_err`` and ``pooled_plugin_std_err``.
"""
names = list(self.variance)
errors = self.plugin_std_error
status = self.inference
payload: dict[str, Any] = {"estimand": names}
if is_labelled(status):
payload["inference"] = [status for _ in names]
payload |= {
"fold_evaluated_psi": [self.fold_evaluated[name].psi for name in names],
"pooled_psi": [self.pooled[name].psi for name in names],
spread_name("cv_std_err", status): [errors[name] for name in names],
spread_name("pooled_std_err", status): [
self.pooled[name].plugin_std_error for name in names
],
"fold_sd": [
float(np.std(self.fold_estimates[name], ddof=1))
if len(self.fold_estimates.get(name, ())) > 1
else float("nan")
for name in names
],
}
return emit_frame(payload, data, backend=self.backend)
def summary(self) -> str:
"""A printable report of the fold-level targeting."""
errors = self.plugin_std_error
status = self.inference
cv_label = spread_name("cv std_err", status)
rows = []
for name, values in self.fold_estimates.items():
spread = f"{np.std(values, ddof=1):.4g}" if len(values) > 1 else "n/a"
rows.append(
[
name,
f"{self.fold_evaluated[name].psi:.5g}",
f"{errors[name]:.4g}",
f"{self.pooled[name].psi:.5g}",
f"{self.pooled[name].plugin_std_error:.4g}",
spread,
f"[{min(values):.5g}, {max(values):.5g}]",
]
)
table = format_table(
[
"estimand",
"fold-evaluated",
cv_label,
"pooled",
spread_name("std_err", status),
"fold sd",
"fold range",
],
rows,
)
header = [
f"Cross-validated targeting over {self.n_folds} folds "
f"(sizes {min(self.fold_sizes)}-{max(self.fold_sizes)})",
"fold-evaluated: evaluated inside each fold and averaged (original CV-TMLE).",
"pooled: this fit's updated validation predictions stitched and evaluated once.",
]
if self.repeats > 1:
header.append(
f"both reports and the {cv_label} combine {self.repeats} draws by their "
"median and split dispersion; the fold columns describe the first."
)
if not supplies_inference(status):
# The fit's own sentence, as ``TMLEResult.summary`` prints it under its table.
header.append(status_record(status).reason)
elif is_labelled(status):
header.append(f"experimental: {status_record(status).shown_reason}")
header.append("")
epsilon_lines = ["", "common fluctuation coefficients:"]
for group, eps in self.epsilon.items():
formatted = "(" + ", ".join(f"{e:.4g}" for e in eps) + ")"
epsilon_lines.append(f" {group}: {formatted}")
if self.fold_epsilon:
epsilon_lines.append("fold-specific fluctuation coefficients:")
for group, per_fold in self.fold_epsilon.items():
formatted = ", ".join(
"(" + ", ".join(f"{e:.4g}" for e in eps) + ")" for eps in per_fold
)
epsilon_lines.append(f" {group}: {formatted}")
return "\n".join([*header, table, *epsilon_lines])
[docs]
@dataclass(frozen=True)
class TMLEResult:
"""Store a fitted point-treatment TMLE result.
Behaves like a mapping from estimand name to
:class:`~cleverly.inference.ParameterEstimate`, and carries the nuisance fits so
that sensitivity and validation analyses can be run without refitting.
Parameters
----------
estimates : dict of str to ParameterEstimate
Estimates keyed by stable parameter alias.
repeats : tuple of RepeatFit
One :class:`~cleverly.estimators._nuisance.RepeatFit` per draw of the
cross-fitting split -- a one-element tuple for an ordinary fit, ``R`` of them
under ``repeats=R``. This is where the nuisances and the fluctuations actually
live; :attr:`nuisance` and :attr:`fluctuations` read through to the first entry,
which is what keeps every analysis written against a single fit working unchanged.
Anything that must account for *all* the draws -- and every analysis that produces
a number, as against one that describes a mechanism, must -- iterates this.
data : CausalData
Validated study data and dataframe backend.
config : TMLEConfig
Normalized implementation configuration.
estimator : Any or None
Fitted estimator that produced this result.
provenance : Provenance or None
Runtime, dependency, and data fingerprints.
simultaneous : SimultaneousBands or None
Joint confidence bands for multi-parameter fits.
bootstrap : BootstrapResult or None
Optional refit-bootstrap draws and failure counts.
intermediate_value : float or None
The level of the intermediate variable this fit targets, or ``None`` for an
ordinary point-treatment fit. It is part of the *estimand*, not a setting: every
estimate here is a controlled direct effect holding ``Z`` at this level, and the
same data yields a different parameter at the other one. See
:mod:`cleverly.estimators.direct_effect`.
extra : dict
Method-specific fitted artifacts.
identified_effect : IdentifiedEffect or None
Causal question and identifying assumptions.
method : EstimationMethod or None
Typed public method configuration.
parameter_keys : dict of str to ParameterKey
Structured identities for reported aliases.
fitted_method : str
Method identity stamped by the estimator; unknown for unstamped artifacts.
solved_corrections : bool
Whether the fit participated in the guarded correction system.
Attributes
----------
assessment_cache : dict
Saved diagnostic and sensitivity outputs. Not a constructor argument, so every
constructed result owns its own cache. A result derived with
:func:`dataclasses.replace` therefore starts empty and cannot serve the original
result's verdicts. Persistence restores the mapping without the constructor.
See Also
--------
cleverly.ParameterEstimate : One entry of :attr:`estimates`.
cleverly.assessment.DiagnosticsFacade : Reached as ``result.diagnostics``.
cleverly.assessment.SensitivityFacade : Reached as ``result.sensitivity``.
cleverly.longitudinal.LongitudinalResult : The same contract for a sequential fit.
Examples
--------
>>> 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=200, seed=0)
>>> 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,
... )
The result is a mapping from estimand alias to
:class:`~cleverly.inference.ParameterEstimate`:
>>> sorted(result.estimates)
['ate']
>>> low, high = result["ate"].ci
>>> low < high
True
The interval bounds are ordered. The default validation battery reads the fitted
artifacts without refitting:
>>> result["ate"].std_error > 0
True
>>> len(result.validate().items) > 0
True
"""
estimates: dict[str, ParameterEstimate]
repeats: tuple[RepeatFit, ...]
data: CausalData
config: TMLEConfig
estimator: Any = None
provenance: Provenance | None = None
simultaneous: SimultaneousBands | None = None
bootstrap: BootstrapResult | None = None
intermediate_value: float | None = None
extra: dict[str, Any] = field(default_factory=dict)
#: Fits made through ``CausalStudy`` record this effect. Direct estimator calls
#: leave it unset.
identified_effect: Any = None
#: The normalized typed method configuration used by ``IdentifiedEffect.estimate``.
method: Any = None
#: Alias-to-structured-key mapping. Routing must read this rather than parse aliases.
parameter_keys: dict[str, Any] = field(default_factory=dict)
fitted_method: str = "unknown"
solved_corrections: bool = False
#: Persistent assessment results, keyed by operation plus normalized arguments.
#: Filling this mapping never changes the fitted parameter or its summary.
#:
#: ``init=False`` because the key records the operation and its arguments and nothing
#: about the result that answered them. ``dataclasses.replace`` used to pass the
#: original's mapping straight through, so a result derived by ``attach_bootstrap`` --
#: which changes ``estimates``, and therefore every sensitivity answer -- served the
#: original's cached verdicts under a different method stamp. A constructed result
#: now owns an empty cache, and joblib persistence restores the saved mapping with the
#: instance state rather than through ``__init__``.
assessment_cache: dict[str, Any] = field(default_factory=dict, init=False)
#: Which family of assessment declarations applies to this result. A class constant
#: rather than a property: the family is a property of the result *type*, and nothing
#: a fit can do changes it. ``cleverly.assessment`` reads it to pick the capability
#: table. The method identity is :attr:`fitted_method`, which is the only name for
#: it, because two names for one value drift.
assessment_family: ClassVar[str] = "point"
# --------------------------------------------------------------- repeats
@property
def nuisance(self) -> NuisanceEstimates:
"""The first draw's nuisance predictions.
The attribute every sensitivity analysis and validation diagnostic was written
against, kept pointing at one object so that a fit with no repeats behaves exactly
as it always has. On a repeated fit this is draw zero of ``R`` and *nothing
else*: it is the right thing to describe a fitted mechanism with, and the wrong
thing to compute a reported number from. Use :attr:`repeats` for the latter.
"""
return self.repeats[0].nuisance
@property
def fluctuations(self) -> dict[str, Fluctuation]:
"""The first draw's solved fluctuations, keyed by target group.
Reads through to :attr:`repeats` for the reason :attr:`nuisance` does, and carries
the same warning: on a repeated fit this is one draw's ``epsilon`` and one draw's
targeted ``Qbar``.
"""
return self.repeats[0].fluctuations
@property
def nuisances(self) -> tuple[NuisanceEstimates, ...]:
"""Every draw's nuisance predictions, in fit order."""
return tuple(repeat.nuisance for repeat in self.repeats)
@property
def n_repeats(self) -> int:
"""How many draws of the cross-fitting split this fit combined."""
return len(self.repeats)
@property
def split_plan(self) -> SplitPlan:
"""Return immutable assignments for every fitted repeat.
Bound to the rows that produced it, through the same data fingerprint
:attr:`provenance` records. A fold label is a *position*, so a plan is only
meaningful for the rows it was realised on: reuse it on reordered rows, or on
rows a refutation replaced, and every label points at a different unit while the
row count still agrees. :meth:`~cleverly.SplitPlan.validate` refuses that, and a
caller who means to reuse the labels on other rows says so with
:meth:`~cleverly.SplitPlan.unbound`.
The plan carries each repeat's :class:`~cleverly.learners.FoldOrigin` when
:func:`~cleverly.learners.random_partition` drew every repeat, and a fit accepts
only such a plan. A stratified or in-sample fit returns a plan whose
:attr:`~cleverly.SplitPlan.provenance` is ``None``: it records the split the fit
ran, and no fit accepts it back.
"""
return SplitPlan.from_folds(
(repeat.folds for repeat in self.repeats),
source_fingerprint=(
None if self.provenance is None else self.provenance.data_fingerprint
),
)
[docs]
def repeat_spread(self) -> dict[str, float]:
r"""Split spread of the estimate across the cross-fitting draws, per estimand.
A *diagnostic*, and emphatically not a standard error. It measures how much the
arbitrary fold assignment moved the answer: the ``R`` draws differ in nothing but
the split, so this is the size of the fold noise a single fit carries silently.
Reported **on the inference scale**, which is the scale
:attr:`~cleverly.ParameterEstimate.std_error` is on: the original outcome scale
for a level or a difference, and the log scale for a ratio. A ``ratio`` estimand
therefore contributes :math:`\mathrm{sd}(\log\hat\psi_r)`. Read against
:attr:`~cleverly.ParameterEstimate.std_error`: a spread that is an appreciable
fraction of the standard error means the split mattered, and one near zero means
the nuisance fits were stable enough that repeating bought little. That is a way
to read the number and not a threshold on it. No rule in this package assigns the
ratio a status, and
:class:`~cleverly.validation.RepeatSpreadRow` reports it without one. A spread taken
on ``psi`` itself would not be comparable with the standard error of an ``rr`` or
an ``or``, and would be wrong by roughly the estimate's own magnitude.
A draw whose ``psi`` is not finite, or whose ratio-scale ``psi`` is not positive,
yields ``nan`` for that estimand rather than ``-inf`` or an exception. There is no
spread on the inference scale to report there, and the caller needs the other
estimands.
What it must not be used for is inference. It says nothing about the *sampling*
variability of the estimand, so it is not an alternative to the reported standard
error. The median variance rule already includes between-draw displacement.
Raises when there is only one draw, since the standard deviation of one number is
not a diagnostic but an artefact.
Returns
-------
dict of str to float
Standard deviation across the cross-fitting draws, per estimand, on the
inference scale of that estimand.
"""
if self.n_repeats < 2:
raise ValueError(
"repeat_spread() describes how much psi moved between draws of the "
f"cross-fitting split, and this fit has {self.n_repeats}. Fit with "
"repeats=2 or more."
)
shared = [name for name in self.estimates if all(name in r.psi for r in self.repeats)]
spreads: dict[str, float] = {}
for name in shared:
values = np.asarray([repeat.psi[name] for repeat in self.repeats], dtype=float)
# Outside the ratio branch, because a draw that is not finite has no spread on
# any scale. Leaving it inside sent an infinite difference-scale draw into
# ``np.std``, which warns and therefore raises under this project's
# ``error::RuntimeWarning`` filter rather than yielding the promised ``nan``.
if not np.all(np.isfinite(values)):
spreads[name] = float("nan")
continue
if self.estimates[name].scale == "ratio":
if np.any(values <= 0.0):
spreads[name] = float("nan")
continue
values = np.log(values)
spreads[name] = float(np.std(values, ddof=1))
return spreads
# ------------------------------------------------------------- accessors
def __getitem__(self, name: str) -> ParameterEstimate:
try:
return self.estimates[name]
except KeyError:
raise KeyError(
f"estimand {name!r} was not requested; available: {list(self.estimates)}"
) from None
def __contains__(self, name: object) -> bool:
return name in self.estimates
def __iter__(self) -> Iterator[str]:
return iter(self.estimates)
@property
def estimate(self) -> ParameterEstimate:
"""The sole parameter estimate, with an explicit refusal for multi-parameter fits."""
return sole_estimate(self.estimates)
@property
def inference_status(self) -> InferenceStatus:
"""The inference status every reported estimate declares.
One status per fit: the estimator stamps it on every estimate it produces, and
:func:`~cleverly.inference.results.inference_status` refuses a mix with
:class:`ValueError`. ``"influence_curve"`` when the fit reports no estimate,
because such a fit refuses nothing.
"""
return reported_status(self.estimates)
[docs]
def psi(self, name: str | None = None) -> float:
"""Return one point estimate.
Parameters
----------
name : str or None
Parameter alias. ``None`` requires a single-parameter result.
Returns
-------
float
Point estimate on its reported scale.
"""
return self.estimate.psi if name is None else self[name].psi
@property
def ate(self) -> ParameterEstimate:
"""The average treatment effect, the most commonly wanted estimand."""
return self["ate"]
@property
def n(self) -> int:
"""Return the number of observations."""
return self.data.n
@property
def influence_curves(self) -> dict[str, FloatArray]:
"""Return influence curves by parameter alias."""
return estimate_curves(self.estimates)
@property
def cv_targeting(self) -> CVTargeting | None:
"""Fold-level detail for fold evaluation or fold-specific targeting.
Carries both the stacked and original fold-evaluated reports whichever one
``result[name]`` was configured to show, so the two can always be compared.
"""
value = self.extra.get("cv_tmle")
return value if isinstance(value, CVTargeting) else None
@property
def ctmle_selection(
self,
) -> CTMLESelection | CTMLEOutcomeAdaptiveFit | CTMLEDerivativeFit | None:
"""The collaborative selection artifact a C-TMLE fit retained, if any.
The sibling of :attr:`cv_targeting` for the ``ctmle`` key, and typed for the same
reason: a reader who reaches into ``extra`` gets ``Any`` and has to re-assert the
type at every call site. ``None`` covers both an ordinary fit and a collaborative
result whose artifact is absent, which the reports distinguish by
:attr:`fitted_method`.
The two classes are imported inside the property because
:mod:`cleverly.estimators.ctmle` imports this module, so the dependency only runs
in the direction the package already has.
"""
from .ctmle import CTMLEDerivativeFit, CTMLEOutcomeAdaptiveFit, CTMLESelection
value = self.extra.get("ctmle")
kinds = CTMLESelection | CTMLEOutcomeAdaptiveFit | CTMLEDerivativeFit
return value if isinstance(value, kinds) else None
# ------------------------------------------------------------- contrasts
[docs]
def covariance(self, names: Sequence[str] | None = None) -> FloatArray:
"""Return joint covariance for selected estimates.
The estimands are *not* independent -- they are functionals of one targeted
distribution and share most of their influence curve -- so a contrast built
from two of them needs this rather than the two variances. Computed from the
influence curves at the right independent unit, so a clustered fit gets the
cluster-level covariance.
Each estimate declares its ``covariance_rule``, and every selection size applies
that one rule. The ordinary rule is the sample covariance of the
curves. On a fold-evaluated CV-TMLE fit, the diagonal therefore differs from the
stored cross-validated variance. The stacked cross-fitted natural-course mean
declares the raw second moment, which is also its stored variance. A selection that
mixes rules, or a raw second-moment selection on a clustered result, raises
``ValueError``.
On a fit whose :attr:`inference_status` supplies no inference, the matrix is the
plug-in covariance of the reported curves. It is a diagnostic and it licenses no
interval, for the reason that status records.
Parameters
----------
names : sequence of str or None
Aliases in output order. ``None`` selects all estimates.
Returns
-------
ndarray
Covariance matrix in the requested order.
"""
refuse_after_repeats(
self.n_repeats,
operation="covariance()",
text=(
"The marginal variance includes between-draw displacement, but no "
"established joint aggregation rule supplies its cross-estimand "
"covariance."
),
reason=RefusalReason.INFERENCE_GAP,
alternative="Fit one split or inspect each entry in result.repeats.",
)
return estimate_covariance(self.estimates, names, cluster=self.data.cluster)
[docs]
def contrast(
self,
function: Callable[[FloatArray], float],
names: Sequence[str],
*,
name: str | None = None,
scale: Scale = "difference",
gradient: Callable[[FloatArray], FloatArray] | None = None,
transform: Transform | None = None,
) -> ParameterEstimate:
r"""A smooth function of several estimands, with correct inference.
Applies the delta method to the *joint* influence curve, so the correlation
between the estimands is handled rather than ignored:
:math:`D_\phi = \nabla\phi(\hat\psi)^\top D`. For example,
``result.contrast(lambda p: p[0] - p[1], ["ey1", "ey0"])`` compares two
means already stored on ``result``.
Pass ``gradient`` when the function's derivative is known in closed form. The
default is a central difference, which is accurate to about ``1e-10`` relative
-- fine for reporting, but not for a test asserting agreement at ``1e-12``.
The result is an ordinary :class:`~cleverly.inference.ParameterEstimate`, so it
carries its own influence curve and can itself be fed back into a contrast.
It inherits the inference status of the estimates it reads, so a contrast of
non-inferential estimates refuses an interval exactly as they do.
``scale="ratio"`` reports a positive contrast with its interval on the log scale.
``transform=`` computes the interval and the test on the scale of a declared
monotone map and maps the interval back with the inverse, as
``ci(contrast = list(f, f_inv, h, fh_grad))`` in R ``drtmle`` 1.1.2 does.
Parameters
----------
function : callable
Smooth scalar function of the selected estimates.
names : sequence of str
Aliases passed to ``function`` in order.
name : str or None
Alias for the derived estimate.
scale : str
Reported scale of the derived estimate. With ``transform=`` it must be
``"level"`` or ``"difference"``.
gradient : callable or None
Analytic gradient. ``None`` uses central differences.
transform : Transform or None
The map the interval and the test are computed on. ``None`` computes them on
the reported scale, or on the log scale for ``scale="ratio"``.
Returns
-------
ParameterEstimate
Derived estimate with influence-curve inference.
"""
refuse_after_repeats(
self.n_repeats,
operation="contrast()",
text=(
"A coordinatewise median does not preserve algebraic identities across estimands."
),
reason=RefusalReason.API_REDIRECTION,
alternative=("Request the contrast as an estimand before fitting, or fit one split."),
)
derived = smooth_contrast(
self.estimates,
function,
names,
n=self.n,
cluster=self.data.cluster,
alpha=self.config.alpha_sig,
name=name,
scale=scale,
gradient=gradient,
transform=transform,
)
return with_derived_bootstrap(
derived,
self.estimates,
self.bootstrap,
names,
function,
cluster=self.data.cluster,
)
[docs]
def ratio(
self,
numerator: str,
denominator: str,
*,
kind: RatioKind = "rr",
name: str | None = None,
) -> ParameterEstimate:
"""The risk ratio or odds ratio of two level estimates of this fit.
The interval is computed on the log scale and exponentiated, as
``lmtp_contrast(type = "rr")`` and ``type = "or"`` in R ``lmtp`` 1.5.4 compute
it. The p-value tests a ratio of one. :meth:`ParameterEstimate.wald_test` tests
any positive null. The ratio inherits the inference status of its inputs.
Parameters
----------
numerator : str
The name of the level in the numerator, such as ``"ey1"``.
denominator : str
The name of the level in the denominator, such as ``"ey0"``.
kind : {"rr", "or"}, default="rr"
A ratio of the levels, or a ratio of their odds. ``"or"`` needs a binary
outcome.
name : str or None, default=None
The name of the derived estimate. ``None`` gives the registered name when
the inputs are an arm and the reference arm, such as ``"rr"`` for
``ratio("ey[1]", "ey[0]")`` on a two-armed fit. Otherwise it gives
``"rr[<numerator> vs <denominator>]"``.
Returns
-------
ParameterEstimate
The ratio, on the ``"ratio"`` scale.
Raises
------
CapabilityError
When ``kind="or"`` and the outcome is not binary.
ValueError
When an input is not a level, or a level is outside the domain of the log.
See Also
--------
cleverly.inference.ParameterEstimate.wald_test : Test the ratio at a null.
Examples
--------
>>> from sklearn.linear_model import LogisticRegression
>>> from cleverly import CausalStudy, CounterfactualMean, PointTreatment
>>> from cleverly.datasets import make_binary_outcome
>>> frame, _ = make_binary_outcome(n=200, seed=0)
>>> study = CausalStudy(
... frame,
... design=PointTreatment(outcome="Y", treatment="A", adjustment=("W1", "W2")),
... )
>>> result = study.identify(CounterfactualMean()).estimate(
... outcome_learner=LogisticRegression(max_iter=1000),
... treatment_learner=LogisticRegression(max_iter=1000),
... cross_fit=False,
... random_state=0,
... )
>>> rr = result.ratio("ey[1]", "ey[0]")
>>> rr.name, rr.scale, rr.psi > 0
('rr', 'ratio', True)
"""
refuse_after_repeats(
self.n_repeats,
operation="ratio()",
text=(
"A coordinatewise median does not preserve algebraic identities across estimands."
),
reason=RefusalReason.API_REDIRECTION,
alternative=("Request the contrast as an estimand before fitting, or fit one split."),
)
if kind == "or":
refuse_odds_ratio(self.config.family)
derived = ratio_contrast(
self.estimates,
numerator,
denominator,
kind=kind,
n=self.n,
cluster=self.data.cluster,
alpha=self.config.alpha_sig,
name=name or self._ratio_name(numerator, denominator, kind),
)
return with_derived_bootstrap(
derived,
self.estimates,
self.bootstrap,
(numerator, denominator),
ratio_function(kind, complement=False),
cluster=self.data.cluster,
)
def _ratio_name(self, numerator: str, denominator: str, kind: str) -> str:
"""The registered name of an arm ratio, or the generic one.
Built forward from the arm labels, never parsed out of the inputs. The registered
name is used when the inputs are the registered level names of two arms. On a
two-armed fit, where the registered name drops the arms, the denominator must
also be the reference arm.
"""
generic = f"{kind}[{numerator} vs {denominator}]"
data = self.data
if self.config.parameter_axis != "arm" or data.is_continuous_treatment:
return generic
collapse = data.is_binary_treatment
# The ``ey`` target names every level ``ey[<label>]``. A two-armed fit also
# names them ``ey1`` and ``ey0`` by arm code.
level_of = {parameter_name("ey", arm=data.arm_label(code)): code for code in data.arm_codes}
if collapse:
level_of.update({"ey1": 1.0, "ey0": 0.0})
top, bottom = level_of.get(numerator), level_of.get(denominator)
if top is None or bottom is None or top == bottom:
return generic
if collapse and bottom != self.config.reference_arm:
return generic
return arm_alias(
kind,
arm=data.arm_label(top),
versus=data.arm_label(bottom),
collapse=collapse,
)
def _names(self, names: Sequence[str] | None) -> tuple[str, ...]:
return select_estimates(self.estimates, names)
# ----------------------------------------------------------- diagnostics
@cached_property
def sensitivity(self) -> SensitivityFacade:
"""Capability-aware sensitivity analyses for this fitted method."""
from ..assessment import SensitivityFacade
return SensitivityFacade(self)
@cached_property
def diagnostics(self) -> DiagnosticsFacade:
"""Unified support, nuisance, score, and refutation diagnostics."""
from ..assessment import DiagnosticsFacade
return DiagnosticsFacade(self)
[docs]
def validate(self) -> ValidationReport:
"""Run the inexpensive method-appropriate checks without refitting.
Returns
-------
ValidationReport
One item per check the fitted method declares, read off stored artifacts.
See Also
--------
cleverly.ValidationReport : The returned validation statuses.
cleverly.assessment.DiagnosticsFacade.run_all : Run the broader diagnostic set.
Examples
--------
>>> 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,
... )
>>> report = result.validate()
>>> len(report.items) > 0
True
"""
from ..assessment import validate_result
return validate_result(self)
[docs]
def assess(
self,
*,
include_refits: bool = False,
include_retargets: bool = False,
arguments: Mapping[str, Mapping[str, Any]] | None = None,
random_state: int | None = None,
) -> AssessmentReport:
"""Run the applicable post-fit assessment battery.
Parameters
----------
include_refits : bool
Run requested operations that refit nuisance models.
include_retargets : bool
Include moderate retargets; cheap E-value retargets run by default.
arguments : mapping or None
Per-operation keyword arguments.
random_state : int or None
Common seed for stochastic refit operations.
Returns
-------
cleverly.AssessmentReport
Validation, diagnostics, and sensitivity in one report.
See Also
--------
TMLEResult.validate : Run the stored-artifact validation checks.
cleverly.assessment.DiagnosticsFacade.run_all : Run diagnostics alone.
Examples
--------
>>> 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)
>>> result = CausalStudy(
... frame,
... design=PointTreatment(
... outcome="Y", treatment="A", adjustment=("W1", "W2", "W3", "W4")
... ),
... ).identify(ATE()).estimate(
... outcome_learner=LinearRegression(),
... treatment_learner=LogisticRegression(max_iter=1000),
... cross_fit=False,
... random_state=0,
... )
>>> result.assess().validation.passed
True
"""
from ..assessment import assess_result
return assess_result(
self,
include_refits=include_refits,
include_retargets=include_retargets,
arguments=arguments,
random_state=random_state,
)
@property
def replayability(self) -> Replayability:
"""Which post-fit operations this in-memory or loaded result can replay."""
from ..assessment import replayability
return replayability(self)
@cached_property
def score_verdict(self) -> ScoreCheck:
"""This fit's own answer to whether its interval is licensed.
The same object ``result.diagnostics.score_equations()`` returns, at the default
tolerance, held here because :meth:`summary` reports it and because a caller
should not have to know which subsystem owns the question. **Derived, never
stored**: it is recomputed from the fluctuations a result carries, so a fit
reloaded from disk answers with its own records rather than with a flag written
at fit time that nothing could check afterwards. Whole-result persistence retains
those records directly instead of rebuilding a partial result graph.
:meth:`__getstate__` is what makes that true across a save, because the memo this
property writes is otherwise part of the artifact.
Free: it reads cached arrays and refits nothing.
"""
from ..validation.score import score_check
return score_check(self)
# ----------------------------------------------------------- persistence
def __getstate__(self) -> dict[str, Any]:
"""Persist the result without the values it derived.
``save`` pickles the result whole, and every
:func:`~functools.cached_property` here writes its answer into ``__dict__``.
``summary()`` reads :attr:`score_verdict`, so the ordinary ``fit``, ``summary``,
``save`` order would bake that verdict into the file. Derived state is not
persisted: the loaded result computes each memo again when it is first read.
:attr:`sensitivity` and :attr:`diagnostics` go with it. Each facade already
drops its own memoized verdicts, so what a stored facade pins is one derived
object graph rather than a stale answer, and rebuilding it costs one
constructor call. Dropping every memo rather than the three by name, because a
fourth is stale in the artifact written the day it is added.
Returns
-------
dict of str to Any
The instance state, without the entries a ``cached_property`` owns.
"""
return without_memos(type(self), self.__dict__)
# ---------------------------------------------------------------- output
[docs]
def save(self, path: Any) -> Any:
"""Write this result to a trusted joblib artifact; see :func:`cleverly.load`.
Loading joblib data can execute arbitrary Python code. Only load files from a
trusted source in a compatible Python and dependency environment. The artifact
records :data:`cleverly.__version__`. A load by a different version warns with
:class:`~cleverly.exceptions.VersionMismatchWarning` and does not migrate the result.
Parameters
----------
path : path-like
Destination file.
Returns
-------
Path
The destination path, as given.
"""
from .serialize import save as _save
return _save(self, path)
[docs]
def to_frame(self) -> Any:
"""Tidy results, one row per estimand, in the caller's dataframe backend.
An ordinary fit emits ``std_err``, ``ci_lower``, ``ci_upper`` and ``p_value``.
A fit whose :attr:`inference_status` supplies no inference emits an
``inference`` column naming that status, and ``plugin_std_err``,
``plugin_interval_lower`` and ``plugin_interval_upper`` in their place. Those
three are a diagnostic and they are not a confidence interval. Exactly one of
``std_err`` and ``inference`` is present in any frame, so a consumer that has to
tell the two apart has a total test. A fit that divided by a declared known
treatment mechanism adds a ``treatment_mechanism`` column that reads ``"known"``.
Returns
-------
dataframe
One row per estimand, in the backend the data arrived in.
"""
rows = [estimate.to_dict() for estimate in self.estimates.values()]
payload: dict[str, Any] = {}
keys: list[str] = []
for row in rows:
for key in row:
if key not in keys:
keys.append(key)
for key in keys:
payload[key] = [row.get(key, None) for row in rows]
if self.simultaneous is not None:
bands = self.simultaneous.bands
payload["simultaneous_lower"] = [
bands.get(name, (None, None))[0] for name in self.estimates
]
payload["simultaneous_upper"] = [
bands.get(name, (None, None))[1] for name in self.estimates
]
if self.intermediate_value is not None:
payload["intermediate"] = [self.intermediate_value] * len(rows)
if self.config.treatment_mechanism == "known":
# Only on a fit that divided by a declared mechanism, so an ordinary frame
# keeps its columns.
payload["treatment_mechanism"] = ["known"] * len(rows)
return self.data.frame_like(payload)
[docs]
def coefficients(self, scale: str = "link") -> Any:
"""A working model's coefficients, on the link scale or exponentiated.
``scale="link"`` reports what the fit estimated: :math:`\\hat\\beta`, with a Wald
interval on the scale the model is linear on. ``scale="ratio"`` reports
:math:`e^{\\hat\\beta}` with the interval exponentiated from that same scale, which
is how the applied literature reports a log- or logit-link marginal structural
model and the reason those links exist here.
**What the exponential means depends on the link, and the two are not
interchangeable.** Under ``link="log"`` a coefficient exponentiates to a risk (or
rate) ratio; under ``link="logit"`` to an *odds* ratio. Reporting one as the other
is a real error rather than a wording preference, so the ``scale`` column of the
returned frame names which it is, row by row. The **intercept** is a third thing
again: :math:`e^{\\beta_0}` is a baseline mean (log) or a baseline odds (logit),
not a ratio of anything, and its p-value tests :math:`\\beta_0 = 0` rather than any
absence of effect. It is reported rather than dropped -- it is a coefficient of
the declared model like the others -- and labelled ``"baseline"``.
Refused on an identity-link fit, where :math:`\\beta` is a risk difference and
:math:`e^\\beta` is not a quantity; and on a fit with no working model at all,
whose parameters are counterfactual means rather than coefficients.
Nothing is re-estimated: the interval and the p-value come from the influence
curve the fit already reported, read on the scale ``scale`` names. This is a
*view*, in the sense :meth:`cleverly.longitudinal.LTMLEResult.curve` is one.
Parameters
----------
scale : {"link", "ratio"}
Coefficient scale to report.
Returns
-------
dataframe
One row per working-model coefficient in the input backend.
"""
if scale not in ("link", "ratio"):
raise ValueError(f"scale must be 'link' or 'ratio'; got {scale!r}")
msm = self.nuisance.msm
if msm is None:
raise CapabilityError(
"this fit has no working model, so it reports counterfactual means rather "
"than coefficients; result.to_frame() is its report.",
reason=RefusalReason.API_REDIRECTION,
alternative=(
"Declare msm= to project those means onto a model whose coefficients "
"are the parameters."
),
)
if scale == "ratio" and msm.link == "identity":
raise CapabilityError(
"an identity-link working model's coefficients are risk differences, and "
"exp() of a difference is not a quantity anybody reports.",
reason=RefusalReason.DIFFERENT_ESTIMAND,
alternative=(
"Declare link='log' for coefficients that are log risk ratios, or "
"link='logit' for log odds ratios, and this view exponentiates those."
),
)
# Term to reported name, composed *forward* through the same rule that named the
# parameters, never split back out of one: a term may legitimately contain a
# bracket, and a name parsed back would then be filed under a term that does not
# exist rather than failing.
from ..targets import parameter_name
ratio = "risk ratio" if msm.link == "log" else "odds ratio"
rows: list[dict[str, Any]] = []
for column, term in enumerate(msm.terms):
estimate = self.estimates[parameter_name("msm", arm=term)]
if scale == "link":
rows.append(estimate.to_dict())
continue
# scale="ratio" with log_psi=beta *is* the exponentiated view: psi becomes
# exp(beta), the interval is built on the log scale and exponentiated, and the
# null moves from zero to one. No arithmetic of its own -- see
# cleverly.inference.ParameterEstimate.ci.
exponentiated = replace(
estimate, psi=float(np.exp(estimate.psi)), log_psi=estimate.psi, scale="ratio"
)
row = exponentiated.to_dict()
# An intercept is found by its column being constant one at every arm, not by
# its name: the term names are the user's own, and "(intercept)" is a
# convention rather than a promise.
row["scale"] = "baseline" if _is_intercept(msm.design[:, :, column]) else ratio
rows.append(row)
return self.data.frame_like({key: [row[key] for row in rows] for key in rows[0]})
[docs]
def influence_frame(self) -> Any:
"""One column per estimand of per-observation influence-curve values.
On a fit whose :attr:`inference_status` supplies no inference, no result shows
the column is the estimator's influence curve, so this frame emits curves and
not inference.
Returns
-------
dataframe
One column per estimand of per-observation influence-curve values.
"""
return self.data.frame_like(
{name: estimate.influence_curve for name, estimate in self.estimates.items()}
)
[docs]
def summary(self, *, protocol: ProtocolDetail = "full") -> str:
"""A printable report of the fit.
Parameters
----------
protocol : {"full", "fingerprint"}
``"full"`` prints every field of the study protocol record. ``"fingerprint"``
prints only its schema and fingerprint, for a reader who has the record already.
Returns
-------
str
A printable report: the estimates, the settings, then the diagnostics.
"""
data = self.data
header = ["Targeted maximum likelihood estimation", "=" * 38]
facts = [
f"n = {data.n}"
+ (
f" ({int(data.observed.sum())} with an observed outcome)"
if data.has_missing_outcome
else ""
)
+ (
f" ({int(data.treatment_recorded.sum())} with a recorded treatment)"
if data.has_missing_treatment
else ""
)
+ f"; covariates = {data.n_covariates}; {_arm_shares(data)}",
]
if self.identified_effect is not None:
facts.extend(self.identified_effect.summary_lines(protocol=protocol))
else:
if data.treatment_adjustment_sources is not None:
# A fit with no study carries no identification record, so the data's own
# declaration prints the same restriction a design records.
facts.append(
"model restrictions: "
+ treatment_mechanism_restriction(data.treatment_adjustment_sources)
)
facts.extend(protocol_summary_lines(None, protocol))
if data.cluster is not None:
facts.append(f"clusters = {_cluster_fact(data)} (cluster-robust variance)")
if data.is_weighted:
report = data.weight_report()
facts.append(
f"observation weights ({report.name or 'weights'}, "
+ ("estimated" if report.estimated else "fixed")
+ f"): effective n = {report.effective_n:.1f}, "
f"design effect = {report.design_effect:.2f}"
)
facts.append(
"estimand: the parameter in the weight-tilted population dP_w = w dP / E[w]; "
"see result.data.weight_report()"
)
if self.intermediate_value is not None:
facts.append(describe_direct_effect(self.intermediate_value, data.intermediate_name))
facts.extend(self.config.describe())
# The selector a collaborative fit stopped at, named in the facts block rather than
# left to the nuisance report. ``describe()`` exists on both artifact classes so a
# caller renders either without discriminating between them.
selection = self.ctmle_selection
if selection is not None:
facts.append(selection.describe())
# Read the same way: ``base`` cannot import ``drtmle``, so the record describes itself.
reduced = self.extra.get("drtmle")
if reduced is not None:
facts.append(reduced.describe())
if data.treatment_adjustment_sources is not None:
facts.append(drtmle_restriction_line(data.treatment_adjustment_sources))
# The same for the learned-rule record, which states the data-adaptive target.
learned = self.extra.get("learned_rule")
if learned is not None:
facts.append(learned.describe())
if self.provenance is not None:
facts.extend(self.provenance.describe())
level = f"{(1 - self.config.alpha_sig) * 100:g}%"
# One fit's estimates carry one status: ``median_estimates`` and
# ``inference_status`` both refuse a mix, so the table is never half refused.
status = self.inference_status
diagnostic = not supplies_inference(status)
t_reference = has_t_reference(self.estimates)
rows = []
record = status_record(status) if diagnostic else None
support_cells = support_column(self.estimates)
if record is not None:
for name, estimate in self.estimates.items():
rows.append([name, f"{estimate.psi:.5g}", f"{estimate.plugin_std_error:.4g}"])
headers = ["estimand", "psi", record.summary_label]
if support_cells is not None:
headers.append("support")
rows = [[*row, cell] for row, cell in zip(rows, support_cells, strict=True)]
table = format_table(headers, rows)
else:
for name, estimate in self.estimates.items():
if not estimate.supplies_inference:
# A parameter status beside an inferential fit status: the row prints
# its point and withholds the rest, and a note below names the status.
rows.append(
[
name,
f"{estimate.psi:.5g}",
"-",
"-",
"-",
*(["-"] if t_reference else []),
]
)
continue
low, high = estimate.ci
rows.append(
[
name,
f"{estimate.psi:.5g}",
f"{estimate.std_error:.4g}",
f"[{low:.5g}, {high:.5g}]",
format_pvalue(estimate.pvalue),
*([reference_label(estimate)] if t_reference else []),
]
)
headers = [
"estimand",
"psi",
"std_err",
f"{level} CI",
"p_value",
*(["df"] if t_reference else []),
]
if support_cells is not None:
headers.append("support")
rows = [[*row, cell] for row, cell in zip(rows, support_cells, strict=True)]
table = format_table(headers, rows)
parts = [*header, *facts, "", table]
notes = support_notes(self.estimates)
if notes:
parts.append("")
parts.extend(notes)
if t_reference:
parts.append("")
assert data.cluster is not None
parts.append(
T_REFERENCE_NOTE.format(
clusters=positive_mass_clusters(data),
fewest=fewest_clusters(
data.cluster,
data.strata,
data.weights if data.is_weighted else None,
),
)
)
if record is not None:
# The whole column is refused, so no ``95% CI`` heading is printed with a "-"
# under it. A dash under that heading still tells a reader an interval is the
# thing that belongs there.
parts.append("")
parts.append(record.summary_note())
else:
withheld: dict[str, list[str]] = {}
for name, estimate in self.estimates.items():
if not estimate.supplies_inference:
withheld.setdefault(estimate.inference, []).append(name)
for flagged, names in withheld.items():
withheld_record = status_record(flagged)
spreads = ", ".join(
f"{name} {self.estimates[name].plugin_std_error:.4g}" for name in names
)
alternative = (
""
if withheld_record.alternative == NO_ALTERNATIVE
else f" {withheld_record.alternative}"
)
parts.append("")
parts.append(f"rows withheld ({', '.join(names)}): {flagged}")
parts.append(
f"{withheld_record.summary_label} (a diagnostic, not a standard error): "
+ spreads
)
parts.append(withheld_record.reason + alternative)
if self.n_repeats > 1:
# Printed under the table rather than left to the scope document, because a
# reader who subtracts two rows of that table gets a third number the fit
# never reported. contrast() refuses for this reason; the table cannot.
parts.append("")
parts.append(
"the report above is coordinatewise: each row is the median of its own "
"draws, so algebraic identities among the rows (ate versus ey1 - ey0, "
"say) do not hold."
)
# Beside the standard error rather than in a separate report, because the
# comparison is the whole content of the number: on its own "0.0065" says
# nothing about whether the split mattered.
parts.append("")
parts.append(
f"split noise -- spread across the {self.n_repeats} draws on the inference "
"scale (the log scale for a ratio), a diagnostic and not a standard error:"
)
for name, value in self.repeat_spread().items():
# The plug-in accessor: this share is a diagnostic ratio, and the fit may
# be one whose standard error the package refuses to name as one.
error = self[name].plugin_std_error
# "-" and a named reason, as every other table in this package renders a
# cell it has no value for. A draw this spread could not be taken over
# printed "nan (nan% of std_err)" here.
if not np.isfinite(value):
parts.append(f" {name:<5s} - (no spread on the inference scale)")
continue
# Named as ``to_frame`` names the same number, so a selector-path fit
# does not describe its diagnostic as a standard error here.
label = spread_name("std_err", status)
share = f"{value / error:.0%} of {label}" if error > 0 else f"{label} unavailable"
parts.append(f" {name:<5s} {value:.4g} ({share})")
if self.simultaneous is None and diagnostic and len(self.estimates) > 1:
# Stated rather than left blank: a reader who asked for bands, or who knows
# they are the default, would otherwise have to guess why none are here.
parts.append("")
parts.append(NO_SIMULTANEOUS_BANDS)
elif self.simultaneous is None and t_reference and len(self.estimates) > 1:
parts.append("")
parts.append(NO_T_REFERENCE_BANDS)
if self.simultaneous is not None:
parts.append("")
parts.append(
f"simultaneous {level} bands (multiplier bootstrap, "
f"critical value {self.simultaneous.critical_value:.3f} vs "
f"{self.simultaneous.pointwise_critical_value:.3f} pointwise)"
f"{' [experimental]' if self.simultaneous.support.label == 'experimental' else ''}:"
)
for name, (low, high) in self.simultaneous.bands.items():
parts.append(f" {name:<5s} [{low:.5g}, {high:.5g}]")
outside = {
name: estimate.inference
for name, estimate in self.estimates.items()
if name not in self.simultaneous.bands
}
if outside:
parts.append(
" the band is joint over the rows above only; it leaves out "
+ ", ".join(f"{name} ({status})" for name, status in outside.items())
+ ", which supply no interval"
)
elif (
not diagnostic
and not t_reference
and len(self.estimates) > 1
and bool(getattr(self.estimator, "simultaneous", False))
and any(not estimate.supplies_inference for estimate in self.estimates.values())
):
# The fit asked for a band, and fewer than two of its estimates supply an
# interval, so it built none. Named rather than left blank.
left_out = ", ".join(
f"{name} ({estimate.inference})"
for name, estimate in self.estimates.items()
if not estimate.supplies_inference
)
parts.append("")
parts.append(
"no simultaneous band: fewer than two estimates supply an interval. "
f"Left out: {left_out}"
)
if self.bootstrap is not None:
parts.append("")
# Named by what each replicate does: it refits the whole estimator on the
# resample. docs/roadmap.md X55 and F2 keep "targeted bootstrap" for a
# distinct procedure.
parts.append(
f"full-refit bootstrap ({self.bootstrap.resampling} resampling, "
f"{self.bootstrap.n_requested - self.bootstrap.n_failed} usable replicates):"
)
for name, estimate in self.estimates.items():
if estimate.bootstrap is None:
continue
low, high = estimate.bootstrap.ci
parts.append(f" {name:<5s} {bootstrap_line(estimate)}")
# Last, and only when it failed. An interval whose score equation is unsolved is
# not a wider interval, it is one the theory does not license, and until now the
# only way to find that out was to know that `validation.score_check()` existed --
# documentation standing in for reporting. Silent on a passing fit deliberately:
# every transcript in the README and the guide is a passing fit, and a line there
# would be noise on the common path and easier to stop reading.
verdict = self.score_verdict
# Under a declared stopping rule the largest score residual, in standard errors, with the
# coverage a Wald interval keeps at it. Empty under the shipped rule, so an ordinary fit
# prints what it always printed.
residual = verdict.rule_line()
if residual:
parts.append("")
parts.append(residual)
if not verdict.passed:
parts.append("")
parts.append(verdict.one_line())
return "\n".join(parts)
def __repr__(self) -> str: # pragma: no cover - debugging aid
return self.summary()
@dataclass(frozen=True)
class TMLEResultSet(Mapping["float | None", TMLEResult]):
"""What :meth:`~cleverly.TMLE.fit` returns: one result per parameter it estimated.
A controlled direct effect is defined *per* value of the intermediate: "the
effect of ``A`` on ``Y`` holding ``Z`` fixed at ``z``" is a different quantity for
each ``z``. Fitting with ``intermediate=`` therefore yields one
:class:`TMLEResult` per level, exactly as R's ``tmle()`` returns a ``tmle.list``.
``fit`` returns this **always**, and an ordinary fit -- no intermediate variable, one
parameter -- is the single-entry case, keyed ``None``. It used to return
``TMLEResult | TMLEResultSet``, which pushed an ``isinstance`` check onto every caller
that could not know in advance which it would get; :class:`CoverageStudy
<cleverly.validation.simulation.CoverageStudy>` carried one, and any user code
branching on ``intermediate=`` carried another. ``None`` is the key rather than an
invented level because it is what the estimator already calls the absence of one.
Use :meth:`single` for the ordinary case. Attribute access is deliberately *not*
forwarded to a lone result: ``__getitem__`` is taken by the level key and cannot also
mean an estimand name, and implicit forwarding would be invisible to a type checker.
The two levels of a controlled direct effect do *not* decompose the total effect into
direct and indirect parts, and their difference is an interaction contrast rather than
a mediated effect. :mod:`cleverly.estimators.direct_effect` writes the parameter
down, derives its influence function, and states the assumptions -- in particular the
one that separates a controlled direct effect from an average treatment effect, that
``W`` suffices to deconfound ``Z -> Y`` as well as ``A -> Y``.
"""
results: dict[float | None, TMLEResult]
intermediate_name: str | None = None
@property
def levels(self) -> tuple[float | None, ...]:
"""The keys, ascending, with the no-intermediate key first."""
return tuple(sorted(self.results, key=_level_order))
def single(self) -> TMLEResult:
"""The sole result.
Raises when the set holds more than one, because there is no defensible choice
between two controlled direct effects: they are different parameters, and picking
one silently is how a script ends up reporting the effect at ``Z = 0`` while its
author believes they asked about ``Z = 1``.
"""
if len(self.results) != 1:
raise KeyError(
f"this fit produced {len(self.results)} results, one per level of "
f"{self.intermediate_name or 'the intermediate'} ({list(self.levels)}), so "
"there is no single one to return. A controlled direct effect is a "
"different parameter at each level -- index the level you mean, e.g. "
f"result[{self.levels[0]!r}]."
)
return next(iter(self.results.values()))
def __getitem__(self, value: float | None) -> TMLEResult:
if isinstance(value, str):
# The likeliest mistake, and worth catching by name: ``TMLEResult`` is indexed
# by estimand and this is indexed by level, so a caller who has not noticed
# which of the two they hold reaches for res["ate"] first. Left to itself,
# float() reports 'could not convert string to float', which says nothing
# about what went wrong.
raise KeyError(
f"a result *set* is indexed by intermediate level, not by estimand: "
f"{value!r} is an estimand name. Use .single()[{value!r}] for an ordinary "
f"fit, or result[<level>][{value!r}] for a controlled direct effect."
)
key = None if value is None else float(value)
try:
return self.results[key]
except KeyError:
raise KeyError(
f"no result for {self.intermediate_name or 'intermediate'} = {value}; "
f"available: {list(self.levels)}"
) from None
def __iter__(self) -> Iterator[float | None]:
return iter(self.levels)
def __len__(self) -> int:
return len(self.results)
def to_frame(self) -> Any:
"""Stacked results with an intermediate-value column."""
frames = [self.results[value].to_frame() for value in self.levels]
if len(frames) == 1:
return frames[0]
import narwhals as nw
stacked = nw.concat([nw.from_native(frame, eager_only=True) for frame in frames])
return stacked.to_native()
def summary(self, *, protocol: ProtocolDetail = "full") -> str:
"""The results, headed by level -- or bare, when there is only one.
A single-entry set prints exactly what its result prints. Heading it
``--- None = ... ---`` would put the internal sentinel in front of every ordinary
fit's output.
Parameters
----------
protocol : {"full", "fingerprint"}
Passed to :meth:`TMLEResult.summary` of each result.
Returns
-------
str
Each result's summary, headed by its level of the intermediate variable.
"""
if len(self.results) == 1 and self.levels[0] is None:
return self.single().summary(protocol=protocol)
blocks = []
for value in self.levels:
label = (
"no intermediate" if value is None else f"{self.intermediate_name} = {value:.0f}"
)
blocks.append(f"--- {label} ---")
blocks.append(self.results[value].summary(protocol=protocol))
return "\n\n".join(blocks)
def __repr__(self) -> str: # pragma: no cover - debugging aid
return self.summary()
def _level_order(value: float | None) -> tuple[int, float]:
"""Sort key that tolerates the ``None`` an ordinary fit is keyed by.
``sorted`` on a mix of ``None`` and floats raises, and a set is only ever homogeneous
by construction -- but relying on that would make the ordering an invariant enforced
nowhere.
"""
return (0, 0.0) if value is None else (1, float(value))
def _is_intercept(column: FloatArray) -> bool:
"""Whether a design column is the constant one at every arm and every row."""
return bool(np.all(np.asarray(column, dtype=float) == 1.0))
def positive_mass_clusters(data: CausalData) -> int:
"""The number of clusters with positive weight mass in the whole fit.
Parameters
----------
data : CausalData
Data with declared clusters.
Returns
-------
int
The count the few-cluster rules read for a marginal estimate.
"""
assert data.cluster is not None
return fewest_clusters(data.cluster, weights=data.weights if data.is_weighted else None)
def _cluster_fact(data: CausalData) -> str:
"""The cluster count, with the size and positive-mass facts a reader needs.
A weighted fit whose clusters hold equal rows and unequal weight mass prints the mass
range in place of the size range. The stratum count applies to a fit with baseline
strata. The facts also name unequal sizes within a stratum and clusters that have
zero weight mass. Equal sizes without strata print the count alone.
"""
cluster = data.cluster
assert cluster is not None
counts = cluster_sizes(cluster)
smallest, largest = int(counts.min()), int(counts.max())
fact = f"{counts.size}"
if smallest != largest:
fact += f", sizes {smallest} to {largest}"
elif data.is_weighted and unequal_cluster_sizes(cluster, data.weights):
mass = cluster_weight_mass(cluster, data.weights)
fact += f", weight mass {mass.min():.4g} to {mass.max():.4g}"
fact += positive_mass_clause(cluster, data.weights if data.is_weighted else None)
if data.has_strata:
assert data.strata is not None
masks = (data.strata == level for level in np.unique(data.strata))
parts = [(cluster[inside], data.weights[inside]) for inside in masks]
uneven = [cluster_sizes(part) for part, _ in parts]
uneven = [part for part in uneven if int(part.min()) != int(part.max())]
if uneven:
fact += (
f", within-stratum sizes {min(int(part.min()) for part in uneven)}"
f" to {max(int(part.max()) for part in uneven)}"
)
elif data.is_weighted:
uneven_mass = [
cluster_weight_mass(part, weight)
for part, weight in parts
if unequal_cluster_sizes(part, weight)
]
if uneven_mass:
low_mass = min(float(part.min()) for part in uneven_mass)
high_mass = max(float(part.max()) for part in uneven_mass)
fact += f", within-stratum weight mass {low_mass:.4g} to {high_mass:.4g}"
stratum_count = fewest_clusters(
cluster, data.strata, data.weights if data.is_weighted else None
)
fact += f", fewest in one stratum {stratum_count}"
return fact
def _arm_shares(data: CausalData) -> str:
"""How the treatment is distributed, for the summary header.
``P(A=1)`` for a binary treatment, unchanged. For more arms that number is not
just uninformative but wrong-looking -- it is the mean of the arm *codes*, so a
three-armed fit reported ``P(A=1) = 0.98`` -- so every arm's share is listed under
its own label instead. A continuous treatment has no arms to take the share of, so
it reports the dose's range: an empty ``arm shares:`` reads as a broken table rather
than as an inapplicable question.
"""
if data.treatment_kind == "vector":
columns = data._component_columns()
parts: list[str] = []
for name, levels in zip(data.component_names, data.component_levels, strict=True):
values = columns[name]
if levels:
listed = ", ".join(
f"{level}={float(np.average(values == level, weights=data.weights)):.3g}"
for level in levels
)
parts.append(f"{name}: {listed}")
else:
dose = np.asarray(values, dtype=float)
parts.append(f"{name}: range [{dose.min():.3g}, {dose.max():.3g}]")
return "vector treatment -- " + "; ".join(parts)
if data.is_continuous_treatment:
dose = np.asarray(data.treatment, dtype=float)
mean = float(np.average(dose, weights=data.weights))
return f"dose: mean {mean:.4g}, range [{dose.min():.3g}, {dose.max():.3g}]"
if data.is_binary_treatment:
return f"P(A=1) = {data.treated_fraction:.4g}"
# ``arm_fractions`` takes the share among the rows whose treatment is recorded.
shares = [
f"{data.arm_label(arm)}={share:.3g}"
for arm, share in zip(data.arm_codes, data.arm_fractions, strict=True)
]
return f"arm shares: {', '.join(shares)}"