"""Causal-question-first public workflow and typed scientific contracts."""
from __future__ import annotations
from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence
from dataclasses import dataclass, field, replace
from typing import TYPE_CHECKING, Any, Literal, Protocol, runtime_checkable
import narwhals as nw
from ._support import combine_support
from ._typing import Family
from .data import CausalData
from .data.causal_data import _vector_components
from .data.known_mechanism import (
KNOWN_CTMLE_REFUSAL,
KNOWN_POLICY_CLASSIFIER_REFUSAL,
probability_columns,
)
from .data.treatment_covariates import (
check_treatment_covariates,
known_and_restricted,
treatment_mechanism_restriction,
)
from .data.validate import RANDOMIZED_INTERCEPT
from .estimators import CTMLE, DRTMLE, TMLE, TMLEResult
from .estimators._strata import check_stratified_targets
from .estimators.composite import missing_treatment_design_refusal
from .exceptions import (
CapabilityError,
CleverlyError,
DataError,
MethodConfigurationError,
Refusal,
RefusalReason,
)
from .inference.multiplier import SimultaneousBands, simultaneous_bands
from .interventions import Incremental, IPSISet, Policy, PolicySet, RegimeSet, Stochastic
from .interventions.base import InterventionKind, refuse_mixed_interventions
from .interventions.learned import LearnedRule, refuse_learned_rule_composition
from .longitudinal import (
IPSI,
LTMLE,
SDR,
LongitudinalData,
LongitudinalResult,
ParametricGFormula,
Repeated,
SequentialGComputation,
SequentialIPW,
)
from .longitudinal import estimator as _sequential
from .longitudinal.estimator import (
BASELINE_MSM_REFUSAL,
BASELINE_TILT_REFUSAL,
IPW_MTP_REFUSAL,
)
from .longitudinal.regimen import declares_mtp, declares_policy, declares_tilt, tilt_nodes
from .methods import (
BASELINE_POINT_REFUSAL,
IPSI_POINT_REFUSAL,
LTMLE_TILT_REFUSAL,
SDR_POINT_REFUSAL,
SDR_TILT_REFUSAL,
CollaborativeTMLEMethod,
DRTMLEMethod,
EstimationMethod,
GComputationMethod,
IPSIMethod,
IPWMethod,
MethodAvailability,
ParametricGFormulaMethod,
Runtime,
SDRMethod,
TMLEMethod,
)
from .msm import MSM, MSMSet, refuse_continuous_msm_mechanisms
from .protocol import ProtocolDetail, StudyProtocol, protocol_summary_lines
from .targets import TARGETS
from .targets.base import (
INTERMEDIATE_MECHANISM,
MISSINGNESS_MECHANISM,
OUTCOME_REGRESSION,
TREATMENT_OBSERVATION_MECHANISM,
Identification,
Interpretation,
arm_alias,
parameter_name,
stratum_alias,
)
from .targets.builtin import (
COMPLETE_OUTCOME_PREFIX,
CONDITIONING_EVENT_PREFIX,
CONSISTENCY_PREFIX,
DR_DESIGN_CONDITIONAL_CLAUSE,
DR_MECHANISM_HEAD,
INTERMEDIATE_CAVEAT_PREFIX,
MISSINGNESS_CAVEAT_PREFIX,
NO_CONFOUNDING_PREFIX,
REFERENCE_ONLY_POSITIVITY_PREFIX,
TWO_ARM_POSITIVITY_PREFIX,
)
from .targets.population_intervention import (
NATURAL_COURSE_TARGET,
POPULATION_INTERVENTION_TARGETS,
)
from .utils.frames import as_frame
if TYPE_CHECKING: # pragma: no cover - typing only
from .assessment import AssessmentReport
__all__ = [
"ATC",
"ATE",
"ATT",
"BackdoorMeanContrast",
"CausalResult",
"CausalStudy",
"ControlledDirectEffect",
"CounterfactualMean",
"Estimand",
"ExplicitAdjustmentProvider",
"IdentificationProvider",
"IdentifiedEffect",
"IncrementalEffect",
"IncrementalMean",
"LearnedRuleValue",
"LongPanel",
"LongitudinalTreatment",
"MSMProjection",
"ModifiedTreatmentPolicy",
"ModifiedTreatmentPolicyEffect",
"NaturalCourseMean",
"OddsRatio",
"ParameterKey",
"PointTreatment",
"PopulationAttributableFraction",
"PopulationAttributableRisk",
"RegimeContrast",
"RegimeMean",
"RiskRatio",
"StudyProtocol",
"TimeToEvent",
]
[docs]
@runtime_checkable
class Estimand(Protocol):
"""Define the public contract for a typed causal question.
Parameters
----------
*args, **kwargs
Present because :func:`typing.runtime_checkable` gives a protocol a synthetic
constructor. A protocol is implemented, not instantiated.
Attributes
----------
name : str
Stable name used in parameter aliases and method dispatch.
definition : str
"""
name: str
@property
def definition(self) -> str:
"""Return the causal quantity in readable notation."""
...
[docs]
@runtime_checkable
class CausalResult(Protocol):
"""Define operations shared by fitted causal-result families.
Concrete point-treatment and longitudinal results provide this interface. Use the
protocol when application code accepts either result family.
Parameters
----------
*args, **kwargs
Present because :func:`typing.runtime_checkable` gives a protocol a synthetic
constructor. A protocol is implemented, not instantiated.
Attributes
----------
estimates : mapping of str to ParameterEstimate
Estimates keyed by stable parameter alias.
identified_effect : IdentifiedEffect
Causal question and identifying assumptions used for the fit.
method : EstimationMethod
Method configuration used for estimation.
parameter_keys : mapping of str to ParameterKey
Structured identities for the reported aliases.
provenance : Provenance
Runtime and dependency information for the fit.
assessment_cache : mapping
Saved diagnostic and sensitivity artifacts.
assessment_family : str
Result-owned capability family declaration.
estimate : ParameterEstimate
The sole parameter estimate. A concrete result raises ``ValueError`` here when it
reports more than one parameter, so ``isinstance`` does not presence-check it.
Use :meth:`psi` or item access to name the parameter you want.
Notes
-----
:func:`typing.runtime_checkable` gives ``isinstance`` a presence check over every
member this class body declares. On Python 3.11 that check reads each member with
:func:`hasattr`, which *calls* a property and swallows only ``AttributeError``. A
member whose read can raise anything else therefore turns ``isinstance`` into an
exception. ``estimate`` is such a member, so this class declares it under
:data:`typing.TYPE_CHECKING`: a type checker still sees it, and the runtime check
does not read it. Add a new member here only when reading it on every supported
result cannot raise.
"""
estimates: Mapping[str, Any]
identified_effect: Any
method: Any
parameter_keys: Mapping[str, Any]
provenance: Any
assessment_cache: Mapping[str, Any]
assessment_family: str
if TYPE_CHECKING: # pragma: no cover - kept out of the runtime presence check
@property
def estimate(self) -> Any:
"""Return the primary parameter estimate."""
...
[docs]
def psi(self, name: str | None = None) -> float:
"""Return one point estimate by alias.
Parameters
----------
name : str or None
Parameter alias. Use ``None`` only when the result has one estimate.
Returns
-------
float
Point estimate on the parameter's reported scale.
"""
...
@property
def influence_curves(self) -> Mapping[str, Any]:
"""Return influence curves keyed by parameter alias."""
...
[docs]
def covariance(self, names: Sequence[str] | None = None) -> Any:
"""Estimate covariance for selected parameters.
Parameters
----------
names : sequence of str or None
Parameter aliases in output order. ``None`` selects every estimate.
Returns
-------
ndarray
Covariance matrix in the requested order.
"""
...
[docs]
def contrast(self, function: Any, names: Sequence[str], **kwargs: Any) -> Any:
"""Apply the delta method to a smooth function of estimates.
Parameters
----------
function : callable
Function from a parameter vector to one scalar.
names : sequence of str
Parameter aliases passed to ``function`` in order.
**kwargs : Any
Options forwarded to the concrete result implementation.
Returns
-------
ParameterEstimate
Transformed estimate and influence-curve inference.
"""
...
[docs]
def ratio(self, numerator: str, denominator: str, **kwargs: Any) -> Any:
"""Return the risk ratio or odds ratio of two level estimates.
Parameters
----------
numerator : str
Alias of the level in the numerator.
denominator : str
Alias of the level in the denominator.
**kwargs : Any
Options forwarded to the concrete result implementation.
Returns
-------
ParameterEstimate
The ratio, with its interval on the log scale.
"""
...
[docs]
def summary(self, *, protocol: ProtocolDetail = "full") -> str:
"""Return a printable fit summary.
Parameters
----------
protocol : {"full", "fingerprint"}
How much of the protocol record the summary prints.
Returns
-------
str
The estimates, the settings, and the recorded facts of the fit.
"""
...
[docs]
def to_frame(self) -> Any:
"""Return estimates in the input dataframe backend."""
...
[docs]
def save(self, path: Any) -> Any:
"""Serialize the fitted result.
Parameters
----------
path : path-like
Destination for the serialized result.
Returns
-------
Path
The destination path, as given.
"""
...
@property
def diagnostics(self) -> Any:
"""Return the diagnostics facade for saved fit artifacts."""
...
[docs]
def validate(self) -> Any:
"""Run validation checks available without new arguments."""
...
[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
-------
AssessmentReport
Validation, diagnostics, and sensitivity in one report.
"""
...
@property
def sensitivity(self) -> Any:
"""Return the sensitivity-analysis facade for this result."""
...
@property
def replayability(self) -> Any:
"""Report which analyses this result can replay."""
...
[docs]
@runtime_checkable
class IdentificationProvider(Protocol):
"""Convert a causal estimand into an observed-data functional.
Implement this protocol to add an identification strategy. The provider must return an
:class:`IdentifiedEffect` with explicit assumptions and nuisance requirements.
Parameters
----------
*args, **kwargs
Present because :func:`typing.runtime_checkable` gives a protocol a synthetic
constructor. A protocol is implemented, not instantiated.
Attributes
----------
name : str
"""
@property
def name(self) -> str:
"""Return the provider name used in reports."""
...
[docs]
def identify(self, study: CausalStudy, estimand: Any) -> IdentifiedEffect:
"""Identify one estimand for a study.
Parameters
----------
study : CausalStudy
Validated observed data and design.
estimand : Estimand
Typed causal quantity to identify.
Returns
-------
IdentifiedEffect
Estimand bound to an observed-data functional and assumptions.
"""
...
[docs]
@dataclass(frozen=True)
class PointTreatment:
"""Declare column roles for a point-treatment study.
Parameters
----------
outcome : str
Outcome column.
treatment : str or sequence of str
Discrete treatment or continuous dose column, or the columns of a vector
treatment. A vector of categorical columns is one categorical treatment over the
observed tuples. A vector with a column in ``continuous_treatment`` takes
``ModifiedTreatmentPolicy`` with :class:`~cleverly.interventions.JointPolicy`
declarations only.
adjustment : sequence of str
Baseline adjustment columns. Supply at least one unless ``randomized=True``.
randomized : bool
Whether the treatment was randomized and an empty adjustment set is intentional.
missingness : str or None
Outcome-observation indicator. One means observed.
intermediate : str or None
Binary intermediate column for a controlled direct effect.
weights : str or None
Probability-weight column.
weights_type : {"probability"}
Interpretation of ``weights``.
weights_estimated : bool
Whether the supplied weights were estimated from these data.
cluster : str or None
Independent-cluster identifier used for variance estimation.
strata : sequence of str
Baseline strata columns, each also an adjustment column. They define
stratum-specific parameters, such as ``ate[S=1]``, reported beside the marginal
ones. :meth:`CausalStudy.identify` refuses a stratum that lacks a treatment arm
beside an incremental estimand, and an ``MSMProjection`` whose design is singular
inside a stratum.
treatment_kind : {"discrete", "continuous"}
Treatment support used to select supported estimands, for a treatment of one
column. A vector treatment keeps ``"discrete"`` and names its doses in
``continuous_treatment``.
outcome_family : {"auto", "gaussian", "binomial"}
Outcome family. ``"auto"`` infers it from observed values.
treatment_missingness : str or None
Treatment-observation indicator. One means the treatment is recorded. It declares
a treatment missing at random. The composite indicator identifies the arm means and
their contrasts on in-sample TMLE and DRTMLE fits, and their point estimates on the
in-sample per-arm outcome-adaptive CTMLE. It also identifies regime means and arm
MSMs on in-sample TMLE fits.
treatment_probabilities : mapping of level to str, or None
The known treatment mechanism of a randomized design (a tuple of ``(level, column)``
pairs, the form the design stores, is also accepted): each treatment level mapped to
the column that holds ``P(A = level | W)``. Every method then divides by these
probabilities in place of an estimate, and the columns are not adjustment
covariates. A design names columns, so the array forms of
:meth:`cleverly.TMLE.fit` are refused here. ``CollaborativeTMLEMethod`` is
unavailable on such a design, because there is no mechanism to select.
treatment_adjustment : sequence of str or None
The adjustment columns that the treatment mechanism reads, a subset of
``adjustment``. ``None`` reads all of them. The outcome regression and every
other nuisance keep the full set. The declaration restricts the observed-data
model and does not change the estimand. The interval is claimed only under that
restriction, and the identification record prints it. An empty sequence declares a
marginal treatment mechanism. To exclude instruments, omit them from
``adjustment`` instead. The declaration is refused beside
``treatment_probabilities``, which supplies a known mechanism, and beside an
incremental estimand, which shifts the fitted propensity score.
continuous_treatment : sequence of str
The columns of a vector ``treatment`` that hold a continuous dose.
treatment_support : Support, mapping of str to Support, or None
The declared conditional support of a continuous dose, as
:meth:`cleverly.TMLE.fit` takes it. A vector treatment maps each continuous
component to its support. A modified treatment policy must assign each row a dose
inside it.
dose_atoms : sequence of float, or mapping of str to sequence of float
The doses of a continuous treatment with a recorded point mass, as
:meth:`cleverly.TMLE.fit` takes them. A vector treatment maps each continuous
component to its atoms.
Attributes
----------
treatment_columns : tuple of str
treatment_name : str
support : {"discrete", "continuous", "vector"}
See Also
--------
LongitudinalTreatment : The same declaration for time-varying treatment.
CausalStudy : What a design is handed to along with the data.
cleverly.datasets.make_linear_ate : A frame with the columns this example names.
Examples
--------
>>> from cleverly import PointTreatment
>>> design = PointTreatment(
... outcome="Y", treatment="A", adjustment=("age", "baseline_score")
... )
>>> design.treatment
'A'
A randomized trial declares that its empty adjustment set is deliberate:
>>> PointTreatment(outcome="Y", treatment="A", randomized=True).adjustment
()
A trial that knows its allocation probabilities names their columns:
>>> PointTreatment(
... outcome="Y",
... treatment="A",
... adjustment=("age",),
... treatment_probabilities={0: "p_control", 1: "p_treated"},
... ).treatment_probabilities
((0, 'p_control'), (1, 'p_treated'))
A design can let the treatment mechanism read a subset of the adjustment set:
>>> PointTreatment(
... outcome="Y",
... treatment="A",
... adjustment=("age", "baseline_score"),
... treatment_adjustment=("age",),
... ).treatment_adjustment
('age',)
"""
outcome: str
treatment: str | tuple[str, ...]
adjustment: Sequence[str] = field(default_factory=tuple)
randomized: bool = False
missingness: str | None = None
intermediate: str | None = None
weights: str | None = None
weights_type: Literal["probability"] = "probability"
weights_estimated: bool = False
cluster: str | None = None
strata: Sequence[str] = field(default_factory=tuple)
treatment_kind: Literal["discrete", "continuous"] = "discrete"
outcome_family: Family = "auto"
treatment_missingness: str | None = None
treatment_probabilities: Any = None
treatment_adjustment: Sequence[str] | None = None
continuous_treatment: Sequence[str] = field(default_factory=tuple)
treatment_support: Any = None
dose_atoms: Sequence[float] | Mapping[str, Sequence[float]] = field(default_factory=tuple)
def __post_init__(self) -> None:
object.__setattr__(self, "adjustment", tuple(self.adjustment))
object.__setattr__(self, "strata", tuple(self.strata))
object.__setattr__(self, "dose_atoms", _point_atoms(self.dose_atoms))
continuous = self.continuous_treatment
object.__setattr__(
self,
"continuous_treatment",
(str(continuous),) if isinstance(continuous, str) else tuple(map(str, continuous)),
)
vector = not isinstance(self.treatment, str)
if vector:
object.__setattr__(self, "treatment", tuple(str(name) for name in self.treatment))
# The container's own rules, so a design and the data it prepares refuse alike.
_vector_components(
self.treatment,
self.continuous_treatment,
None if vector and self.treatment_kind == "discrete" else self.treatment_kind,
)
if self.treatment_probabilities is not None:
columns = probability_columns(_stored_pairs(self.treatment_probabilities))
if columns is None:
raise DataError(
"PointTreatment(treatment_probabilities=) maps each treatment level to "
"the column that holds its probability, as every other design field "
"names columns. Put the probabilities in the frame, or pass the arrays "
"to TMLE.fit(treatment_probabilities=)."
)
object.__setattr__(self, "treatment_probabilities", columns)
if not self.outcome or not self.treatment:
raise DataError("outcome and treatment must be non-empty column names")
if self.outcome in self.treatment_columns:
raise DataError("outcome and treatment must name different columns")
if len(set(self.adjustment)) != len(self.adjustment):
raise DataError(f"adjustment contains duplicate columns: {list(self.adjustment)}")
if not self.adjustment and not self.randomized:
raise DataError(
"declare a non-empty adjustment set or randomized=True; an empty set is "
"an identification claim, not a nuisance-model default"
)
if self.treatment_adjustment is not None:
object.__setattr__(
self,
"treatment_adjustment",
check_treatment_covariates(
self.treatment_adjustment,
self.adjustment,
argument="treatment_adjustment",
base="adjustment",
),
)
if self.treatment_probabilities is not None:
raise known_and_restricted("treatment_adjustment")
@property
def treatment_columns(self) -> tuple[str, ...]:
"""The treatment's columns: one, or the components of a vector."""
return (self.treatment,) if isinstance(self.treatment, str) else self.treatment
@property
def treatment_name(self) -> str:
"""The treatment's name: its column, or its components joined by ``+``."""
return "+".join(self.treatment_columns)
@property
def support(self) -> Literal["discrete", "continuous", "vector"]:
"""What the prepared treatment is, as :attr:`CausalData.treatment_kind` records it.
``"vector"`` is a vector with a continuous component, which has no arms. An
all-categorical vector is ``"discrete"``: one categorical treatment over its tuples.
"""
if isinstance(self.treatment, str):
return self.treatment_kind
return "vector" if self.continuous_treatment else "discrete"
[docs]
def prepare(self, data: Any) -> CausalData:
"""Validate data and encode the declared column roles.
Parameters
----------
data : dataframe or CausalData
Pandas, Polars, or prepared causal data.
Returns
-------
CausalData
Validated internal representation.
"""
if isinstance(data, CausalData):
self._check_prepared(data)
return data
covariates = self.adjustment
if self.randomized and not covariates:
frame = as_frame(data)
if RANDOMIZED_INTERCEPT in frame.columns:
raise DataError(
f"reserved internal column {RANDOMIZED_INTERCEPT!r} is already present; "
"rename it"
)
data = frame.with_columns(nw.lit(0.0).alias(RANDOMIZED_INTERCEPT)).to_native()
covariates = (RANDOMIZED_INTERCEPT,)
vector = not isinstance(self.treatment, str)
return CausalData.from_frame(
data,
outcome=self.outcome,
treatment=list(self.treatment) if vector else self.treatment,
covariates=covariates,
delta=self.missingness,
intermediate=self.intermediate,
weights=self.weights,
weights_type=self.weights_type,
weights_estimated=self.weights_estimated,
id=self.cluster,
strata=self.strata,
family=self.outcome_family,
treatment_kind=None if vector else self.treatment_kind,
continuous_treatment=self.continuous_treatment,
treatment_delta=self.treatment_missingness,
treatment_probabilities=(
None if self.treatment_probabilities is None else dict(self.treatment_probabilities)
),
treatment_covariates=self.treatment_adjustment,
).with_dose_declarations(self.treatment_support, _declared_atoms(self.dose_atoms))
def _check_prepared(self, data: CausalData) -> None:
"""Reconcile an already-built container with the roles this design declares.
A ``CausalData`` handed straight to ``CausalStudy`` used to be adopted as-is, so
every field here was a claim nothing checked. The design is what
``IdentifiedEffect.functional`` records and what ``summary()`` prints, so a design
naming ``adjustment=("Z1",)`` over a container fitted on ``W1..W3`` reported an
adjustment set no estimate came from -- and, since this PR, persisted it.
"""
expected: tuple[tuple[str, Any, Any], ...] = (
("outcome", self.outcome, data.outcome_name),
("treatment", self.treatment_name, data.treatment_name),
(
"continuous_treatment",
tuple(sorted(self.continuous_treatment)),
tuple(
sorted(
name
for name, levels in zip(
data.component_names, data.component_levels, strict=True
)
if not levels
)
),
),
("missingness", self.missingness, data.delta_name),
("treatment_missingness", self.treatment_missingness, data.treatment_delta_name),
("intermediate", self.intermediate, data.intermediate_name),
("cluster", self.cluster, data.cluster_name),
("weights", self.weights, data.weights_name),
("strata", tuple(self.strata), tuple(data.strata_names)),
("treatment_kind", self.support, data.treatment_kind),
(
"treatment_support",
_point_support(self.treatment_support),
_point_support(data.treatment_support),
),
("dose_atoms", _point_atoms(self.dose_atoms), _point_atoms(data.dose_atoms)),
# Not role names, but they change what the reported variance means, so a
# disagreement here is as consequential as a mis-named column.
("weights_type", self.weights_type, data.weight_spec.kind),
("weights_estimated", self.weights_estimated, data.weight_spec.estimated),
)
for role, declared, held in expected:
if declared != held:
raise DataError(
f"the supplied CausalData was built with {role}={held!r}, but this design "
f"declares {role}={declared!r}; build the data from this design, or "
"correct the design"
)
held = None if data.known_treatment is None else dict(data.known_treatment.columns)
declared = (
None if self.treatment_probabilities is None else dict(self.treatment_probabilities)
)
if held != declared:
raise DataError(
f"the supplied CausalData declares treatment_probabilities={held!r}, but this "
f"design declares treatment_probabilities={declared!r}; build the data from "
"this design, or correct the design"
)
# ``outcome_family`` is inferred when the design does not state one, so only a
# stated disagreement is a disagreement.
if self.outcome_family != "auto" and self.outcome_family != data.family:
raise DataError(
f"the supplied CausalData has family={data.family!r}, but this design "
f"declares outcome_family={self.outcome_family!r}"
)
# ``covariate_names`` is post-encoding, so compare the columns the caller named:
# a categorical adjustment variable arrives as several generated columns, and a
# degenerate one may have been dropped entirely.
generated = data.indicator_sources
sources = {generated.get(name, name) for name in data.covariate_names}
sources.update(data.dropped_covariates)
declared_adjustment = set(self.adjustment) or {RANDOMIZED_INTERCEPT}
if sources != declared_adjustment:
raise DataError(
f"the supplied CausalData adjusts for {sorted(sources)}, but this design "
f"declares {sorted(declared_adjustment)}; an identification claim and the "
"data it is claimed about have to be the same set"
)
# The mechanism's covariates are a restriction the interval depends on, so a
# container that reads another subset is refused like any other mismatch. A
# declaration of the whole adjustment set is no restriction, and a covariate that
# was dropped as degenerate is in neither set.
dropped = set(data.dropped_covariates)
restricted: set[str] | None = (
None if self.treatment_adjustment is None else set(self.treatment_adjustment) - dropped
)
if restricted == set(self.adjustment) - dropped:
restricted = None
used = data.treatment_adjustment_sources
if (restricted is None) != (used is None) or (
used is not None and restricted is not None and set(used) != restricted
):
raise DataError(
f"the supplied CausalData has treatment_adjustment={used!r}, but this design "
f"declares treatment_adjustment={self.treatment_adjustment!r}; build the data "
"from this design, or correct the design"
)
def _stored_pairs(supplied: Any) -> Any:
"""Read the stored form of a mapping field, a tuple of key and value pairs, as a mapping.
A frozen design stores each declared mapping as a tuple of pairs, so that the design is
hashable. ``dataclasses.replace`` hands every stored field back to ``__init__``, so each
normalizer has to read the form it stores. Any other value is returned unchanged. Pairs
that repeat a key, or hold an unhashable key, say what a mapping cannot say, so they raise
:class:`~cleverly.exceptions.DataError`.
"""
if isinstance(supplied, tuple) and all(
isinstance(pair, tuple) and len(pair) == 2 for pair in supplied
):
keys = [key for key, _ in supplied]
try:
distinct = len(set(keys))
except TypeError:
raise DataError(
f"a mechanism declared as pairs needs hashable keys; got {keys!r}"
) from None
if distinct != len(keys):
repeated = next(key for key in keys if keys.count(key) > 1)
raise DataError(
f"a mechanism declared as pairs repeats the key {repeated!r}; a mapping "
"holds one column for each key"
)
return dict(supplied)
return supplied
def _node_probability_columns(supplied: Any) -> tuple[tuple[str, tuple[tuple[Any, str], ...]], ...]:
"""A longitudinal design's treatment declaration, as hashable column pairs by node."""
pairs = None
supplied = _stored_pairs(supplied)
if isinstance(supplied, Mapping):
pairs = tuple(
(str(node), probability_columns(_stored_pairs(per_node)))
for node, per_node in supplied.items()
)
if pairs is None or any(columns is None for _, columns in pairs):
raise DataError(
"LongitudinalTreatment(treatment_probabilities=) maps each treatment node to a "
"mapping of its levels to the columns that hold their probabilities, as every "
"other design field names columns. Put the probabilities in the frame, or pass "
"the arrays to LTMLE.fit(treatment_probabilities=)."
)
return tuple((node, columns) for node, columns in pairs if columns is not None)
def _censoring_probability_columns(supplied: Any) -> tuple[tuple[str, str], ...]:
"""A longitudinal design's censoring declaration, as hashable column pairs."""
supplied = _stored_pairs(supplied)
if not isinstance(supplied, Mapping) or not all(
isinstance(column, str) for column in supplied.values()
):
raise DataError(
"LongitudinalTreatment(censoring_probabilities=) maps each censoring column to the "
"column that holds its retention probability. Put the probabilities in the frame, "
"or pass the arrays to LTMLE.fit(censoring_probabilities=)."
)
return tuple((str(node), str(column)) for node, column in supplied.items())
def _held_node_columns(
data: LongitudinalData,
) -> tuple[tuple[str, tuple[tuple[Any, str], ...]], ...] | None:
"""The treatment column declaration a prepared container holds, as the design spells it."""
known = data.known_mechanisms
if known is None or not known.treatment_columns:
return None
return tuple(
(name, pairs)
for name, pairs in zip(data.treatment_names, known.treatment_columns, strict=False)
if pairs
)
def _held_censoring_columns(data: LongitudinalData) -> tuple[tuple[str, str], ...] | None:
"""The censoring column declaration a prepared container holds, as the design spells it."""
known = data.known_mechanisms
if known is None or not known.censoring_columns:
return None
return tuple(
(name, column)
for name, column in zip(data.censoring_names, known.censoring_columns, strict=False)
if column is not None
)
[docs]
@dataclass(frozen=True)
class LongitudinalTreatment:
"""Declare time-ordered roles for a longitudinal study.
Parameters
----------
outcome : str, sequence of str, mapping of str to sequence of str, or Repeated
End-of-study outcome, survival event nodes, competing-risk event nodes, or a
:class:`~cleverly.longitudinal.Repeated` outcome measured at every node.
treatment : str or sequence of str
Treatment nodes in time order. One column name declares one baseline decision
held over every node. The node count then comes from a survival or repeated
outcome, else from ``censoring``, else from ``time_varying``, else it is one.
baseline : sequence of str
Covariates observed before the first treatment. The list names at least one column.
time_varying : sequence of sequence of str or None
Covariate blocks observed before each treatment node.
censoring : sequence of str or None
Observation indicators in time order. One means observed.
cluster : str or None
Independent-cluster identifier.
weights : str or None
Probability-weight column.
weights_type : {"probability"}
Interpretation of ``weights``.
weights_estimated : bool
Whether the supplied weights were estimated from these data.
outcome_family : {"auto", "gaussian", "binomial"}
Outcome family used by the sequential regressions.
continuous_treatment : sequence of str, or mapping of str to Support or None
The treatment nodes that hold a continuous dose. Each such node takes a modified
treatment policy, such as :class:`~cleverly.interventions.Shift`. A mapping
also declares the conditional support of each dose, as
:meth:`cleverly.longitudinal.LTMLE.fit` takes it.
dose_atoms : mapping of str to sequence of float, or None
The doses of each continuous node with a recorded point mass, as
:meth:`cleverly.longitudinal.LTMLE.fit` takes them.
treatment_probabilities : mapping or None
The known treatment mechanism (a tuple of pairs, the form the design stores, is also
accepted) of a sequentially randomized design: each treatment
node mapped to a mapping of its levels to the column of
``P(A_t = level | observed past)``. A declared node fits no treatment learner.
censoring_probabilities : mapping or None
The known retention probabilities (a tuple of pairs is also accepted): each censoring
column mapped to the column of ``P(C_t = 1 | observed past)``. A declared node fits no
censoring learner.
See Also
--------
PointTreatment : The same declaration for treatment given once.
TimeToEvent : A baseline treatment with one time and one event column.
CausalStudy : What a design is handed to along with the data.
cleverly.datasets.make_longitudinal : A frame with the nodes this example names.
Examples
--------
Two treatment nodes, with the covariates measured before each one and an observation
indicator per node:
>>> from cleverly import LongitudinalTreatment
>>> design = LongitudinalTreatment(
... outcome="Y",
... treatment=["A1", "A2"],
... baseline=["W1", "W2"],
... time_varying=[[], ["L2"]],
... censoring=["C1", "C2"],
... )
>>> design.treatment
('A1', 'A2')
``time_varying`` is indexed by treatment node, so the first block is empty here: no
covariate is measured between baseline and the first treatment.
>>> design.time_varying
((), ('L2',))
"""
outcome: str | Sequence[str] | Mapping[str, Sequence[str]] | Repeated
treatment: str | Sequence[str]
baseline: Sequence[str]
time_varying: Sequence[Sequence[str]] | None = None
censoring: Sequence[str] | None = None
cluster: str | None = None
weights: str | None = None
weights_type: Literal["probability"] = "probability"
weights_estimated: bool = False
outcome_family: Family = "auto"
continuous_treatment: Sequence[str] | Mapping[str, Any] = ()
treatment_probabilities: Any = None
censoring_probabilities: Any = None
dose_atoms: Mapping[str, Sequence[float]] | None = None
def __post_init__(self) -> None:
continuous = self.continuous_treatment
object.__setattr__(
self,
"continuous_treatment",
tuple((str(name), support) for name, support in continuous.items())
if isinstance(continuous, Mapping)
else tuple(continuous),
)
if self.dose_atoms is not None:
object.__setattr__(
self,
"dose_atoms",
tuple(
(str(name), tuple(sorted(float(atom) for atom in atoms)))
for name, atoms in dict(self.dose_atoms).items()
),
)
if self.treatment_probabilities is not None:
object.__setattr__(
self,
"treatment_probabilities",
_node_probability_columns(self.treatment_probabilities),
)
if self.censoring_probabilities is not None:
object.__setattr__(
self,
"censoring_probabilities",
_censoring_probability_columns(self.censoring_probabilities),
)
if not isinstance(self.treatment, str):
object.__setattr__(self, "treatment", tuple(self.treatment))
object.__setattr__(self, "baseline", tuple(self.baseline))
if self.time_varying is not None:
object.__setattr__(
self, "time_varying", tuple(tuple(block) for block in self.time_varying)
)
if self.censoring is not None:
object.__setattr__(self, "censoring", tuple(self.censoring))
if not self.treatment:
raise DataError("a longitudinal design needs at least one treatment node")
[docs]
def prepare(self, data: Any) -> LongitudinalData:
"""Validate data and construct its ordered longitudinal representation.
Parameters
----------
data : dataframe or LongitudinalData
Pandas, Polars, or prepared longitudinal data.
Returns
-------
LongitudinalData
Validated time-ordered representation.
"""
if isinstance(data, LongitudinalData):
self._check_prepared(data)
return data
return LongitudinalData.from_frame(
data,
outcome=self.outcome,
treatment=self.treatment,
baseline=self.baseline,
time_varying=self.time_varying,
censoring=self.censoring,
id=self.cluster,
weights=self.weights,
weights_type=self.weights_type,
weights_estimated=self.weights_estimated,
family=self.outcome_family,
continuous_treatment=_declared_continuous(self.continuous_treatment),
treatment_probabilities=(
None
if self.treatment_probabilities is None
else {node: dict(pairs) for node, pairs in self.treatment_probabilities}
),
censoring_probabilities=(
None if self.censoring_probabilities is None else dict(self.censoring_probabilities)
),
dose_atoms=None if self.dose_atoms is None else dict(self.dose_atoms),
)
def _outcome_roles(self) -> tuple[tuple[str, ...], tuple[str, ...]]:
"""The event nodes and cause labels this outcome declaration implies.
One field spells three designs: a name is an end-of-study outcome, a sequence is a
survival process, and a mapping of cause to sequence is competing risks. The
container records them as ``event_names`` plus ``cause_labels``, so composing them
forward is what makes the three comparable without inspecting arrays.
"""
if isinstance(self.outcome, (str, Repeated)):
return (), ()
if isinstance(self.outcome, Mapping):
causes = tuple(self.outcome)
return tuple(node for cause in causes for node in self.outcome[cause]), causes
return tuple(self.outcome), ()
def _check_prepared(self, data: LongitudinalData) -> None:
"""Reconcile an already-built container with the roles this design declares.
The point-treatment reasoning applies unchanged, and the node ordering makes it
sharper: ``functional.adjustment`` records this design's baseline, and a sequential
fit's history is not recoverable from the arrays afterwards.
"""
events, causes = self._outcome_roles()
empty: tuple[str, ...] = ()
held = isinstance(self.treatment, str)
# A vector node is declared as its list of columns and named by their join.
declared_treatment = (
(self.treatment,)
if held
else tuple(node if isinstance(node, str) else "+".join(node) for node in self.treatment)
)
blocks = (
tuple(empty for _ in range(data.n_times if held else len(declared_treatment)))
if self.time_varying is None
else tuple(tuple(block) for block in self.time_varying)
)
repeated = self.outcome.columns if isinstance(self.outcome, Repeated) else ()
# The container records Y_s as the first covariate of node s + 1, unless the covariate
# screen dropped a constant copy. The repeated nodes are compared on their own below,
# so they are removed from both sides of the covariate comparison.
unmatched = {*repeated, *data.repeated_names}
blocks = tuple(tuple(name for name in block if name not in unmatched) for block in blocks)
prepared_blocks = tuple(
tuple(name for name in block if name not in unmatched)
for block in data.time_varying_names
)
expected: tuple[tuple[str, Any, Any], ...] = (
("treatment", declared_treatment, tuple(data.treatment_names)),
("held treatment", held, data.treatment_held),
("time grid", None, data.time_grid),
("baseline", tuple(self.baseline), tuple(data.baseline_names)),
("time_varying", blocks, prepared_blocks),
("censoring", tuple(self.censoring or ()), tuple(data.censoring_names)),
("cluster", self.cluster, data.cluster_name),
("weights", self.weights, data.weights_name),
("weights_type", self.weights_type, data.weight_spec.kind),
("weights_estimated", self.weights_estimated, data.weight_spec.estimated),
("outcome event nodes", events, tuple(data.event_names)),
("repeated outcome nodes", repeated, tuple(data.repeated_names)),
("outcome causes", causes, tuple(data.cause_labels)),
(
"continuous_treatment",
_continuous_names(self.continuous_treatment),
data.continuous_columns,
),
(
"the declared dose supports",
_named_declarations(self.continuous_treatment),
_held_supports(data),
),
(
"dose_atoms",
dict(self.dose_atoms or ()),
_held_atoms(data),
),
(
"treatment_probabilities",
self.treatment_probabilities,
_held_node_columns(data),
),
(
"censoring_probabilities",
self.censoring_probabilities,
_held_censoring_columns(data),
),
)
for role, declared, held in expected:
if declared != held:
raise DataError(
f"the supplied LongitudinalData was built with {role}={held!r}, but this "
f"design declares {role}={declared!r}; build the data from this design, "
"or correct the design"
)
if not events and isinstance(self.outcome, str) and self.outcome != data.outcome_name:
raise DataError(
f"the supplied LongitudinalData was built with outcome={data.outcome_name!r}, "
f"but this design declares outcome={self.outcome!r}"
)
if self.outcome_family != "auto" and self.outcome_family != data.family:
raise DataError(
f"the supplied LongitudinalData has family={data.family!r}, but this design "
f"declares outcome_family={self.outcome_family!r}"
)
[docs]
@dataclass(frozen=True)
class TimeToEvent:
r"""Declare a baseline treatment with one follow-up time and one event code per unit.
The data hold one row per unit. The time is binned onto a grid
:math:`0 = g_0 < g_1 < \dots < g_K`, node :math:`k` covers :math:`(g_{k-1}, g_k]`, and
the treatment is one decision held over every node.
:meth:`cleverly.longitudinal.LongitudinalData.from_time_to_event` states the binning
rules and the refusals. The estimand is a risk or a cumulative incidence at the grid
times, as for a survival :class:`LongitudinalTreatment`. At a grid time before which
no follower of a regimen had the event, the risk is zero and reports no interval
(status ``"constant_node_plugin"``).
Parameters
----------
time : str
The follow-up time column. Every time is a positive finite number.
event : str
The event-code column: ``0`` for censored, and a positive integer per cause.
treatment : str
The baseline treatment column.
baseline : sequence of str
Covariates measured before the treatment. The list names at least one column.
grid : sequence of float or None
The grid times, which are the visit times. ``None`` needs integer times and uses
:math:`g_k = k` up to the largest event time. The fit costs :math:`K(K+1)/2`
sequential regressions per regimen and cause, so a default grid on day-level times
can hold thousands of nodes.
causes : mapping of int to str or None
The label of each nonzero event code. ``None`` reads codes ``{0, 1}`` as one event
and any other codes as competing causes labelled ``str(code)``.
cluster : str or None
Independent-cluster identifier.
weights : str or None
Probability-weight column.
weights_type : {"probability"}
Interpretation of ``weights``.
weights_estimated : bool
Whether the supplied weights were estimated from these data.
time_varying : None
Not accepted. One row per unit holds no covariate measured after baseline, so any
other value raises :class:`~cleverly.exceptions.CapabilityError` with code
``api_redirection``.
continuous_treatment : bool
Whether the treatment is a continuous dose. A regimen is then a modified treatment
policy at baseline.
treatment_probabilities : mapping of level to str, or None
The known baseline treatment mechanism (a tuple of ``(level, column)`` pairs is also
accepted) of a randomized design: each treatment level
mapped to the column that holds :math:`P(A = \text{level} \mid W)`. The fit then
divides by these probabilities and fits no treatment learner. A design names
columns, so the array forms of :meth:`cleverly.longitudinal.LTMLE.fit` are refused
here.
treatment_support : Support or None
The declared conditional support of a continuous treatment.
dose_atoms : sequence of float
The doses of a continuous treatment with a recorded point mass.
Attributes
----------
outcome : str
See Also
--------
LongitudinalTreatment : The wide layout, which also takes post-baseline covariates.
CausalStudy : What a design is handed to along with the data.
cleverly.longitudinal.LongitudinalData.from_time_to_event : The container it builds.
Examples
--------
>>> from cleverly import TimeToEvent
>>> design = TimeToEvent(
... time="ftime", event="ftype", treatment="arm", baseline=["age"], grid=[30, 60, 90]
... )
>>> design.grid
(30.0, 60.0, 90.0)
"""
time: str
event: str
treatment: str
baseline: Sequence[str]
grid: Sequence[float] | None = None
causes: Mapping[int, str] | None = None
cluster: str | None = None
weights: str | None = None
weights_type: Literal["probability"] = "probability"
weights_estimated: bool = False
time_varying: None = None
continuous_treatment: bool = False
treatment_probabilities: Any = None
treatment_support: Any = None
dose_atoms: Sequence[float] = ()
def __post_init__(self) -> None:
object.__setattr__(self, "dose_atoms", tuple(float(atom) for atom in self.dose_atoms))
if self.treatment_probabilities is not None:
columns = probability_columns(_stored_pairs(self.treatment_probabilities))
if columns is None:
raise DataError(
"TimeToEvent(treatment_probabilities=) maps each treatment level to the "
"column that holds its probability, as every other design field names "
"columns. Put the probabilities in the frame, or pass the arrays to "
"LTMLE.fit(treatment_probabilities=)."
)
object.__setattr__(self, "treatment_probabilities", columns)
if self.time_varying is not None:
raise CapabilityError(
"TimeToEvent reads one row per unit, so it has no covariate measured after "
"baseline.",
reason=RefusalReason.API_REDIRECTION,
alternative=(
"Use LongitudinalTreatment(treatment='A', time_varying=...) on the wide layout."
),
)
if not isinstance(self.treatment, str):
raise TypeError(
"TimeToEvent takes one baseline treatment column; a treatment that changes "
"over time needs LongitudinalTreatment on the wide layout"
)
object.__setattr__(self, "baseline", tuple(self.baseline))
if self.grid is not None:
object.__setattr__(self, "grid", tuple(float(g) for g in self.grid))
if self.causes is not None:
object.__setattr__(self, "causes", dict(self.causes))
@property
def outcome(self) -> str:
"""The event-code column, which the identified functional names as the outcome."""
return self.event
[docs]
def prepare(self, data: Any) -> LongitudinalData:
"""Validate data and construct its held survival representation.
Parameters
----------
data : dataframe or LongitudinalData
Pandas, Polars, or prepared longitudinal data.
Returns
-------
LongitudinalData
The held container with its time grid.
"""
if isinstance(data, LongitudinalData):
self._check_prepared(data)
return data
built = LongitudinalData.from_time_to_event(
data,
time=self.time,
event=self.event,
treatment=self.treatment,
baseline=self.baseline,
grid=self.grid,
causes=self.causes,
id=self.cluster,
weights=self.weights,
weights_type=self.weights_type,
weights_estimated=self.weights_estimated,
continuous_treatment=self.continuous_treatment,
treatment_probabilities=(
None if self.treatment_probabilities is None else dict(self.treatment_probabilities)
),
)
return _declare_every_node(built, self.treatment_support, self.dose_atoms)
def _check_prepared(self, data: LongitudinalData) -> None:
"""Reconcile an already-built container with the roles this design declares."""
held_columns = dict(_held_node_columns(data) or ()).get(self.treatment)
expected: tuple[tuple[str, Any, Any], ...] = (
("treatment_probabilities", self.treatment_probabilities, held_columns),
("time", self.time, data.time_name),
("event", self.event, data.event_name),
("treatment", (self.treatment,), tuple(data.treatment_names)),
("held treatment", True, data.treatment_held),
("baseline", tuple(self.baseline), tuple(data.baseline_names)),
("cluster", self.cluster, data.cluster_name),
("weights", self.weights, data.weights_name),
("weights_type", self.weights_type, data.weight_spec.kind),
("weights_estimated", self.weights_estimated, data.weight_spec.estimated),
("continuous_treatment", self.continuous_treatment, data.is_continuous_node(1)),
("treatment_support", self.treatment_support, data.node_support(1)),
("dose_atoms", tuple(sorted(self.dose_atoms)), data.node_atoms(1)),
)
for role, declared, held in expected:
if declared != held:
raise DataError(
f"the supplied LongitudinalData was built with {role}={held!r}, but this "
f"design declares {role}={declared!r}; build the data from this design, "
"or correct the design"
)
grid = self.grid if self.grid is not None else _default_grid(data)
if tuple(grid) != data.time_grid:
declared_grid = "the default grid " if self.grid is None else ""
raise DataError(
f"the supplied LongitudinalData was built with grid={data.time_grid!r}, but "
f"this design declares {declared_grid}grid={tuple(grid)!r}"
)
held_causes = (
dict(zip(data.event_codes, data.cause_labels, strict=True))
if data.cause_labels
else None
)
declared_causes = (
dict(self.causes)
if self.causes is not None
else None
if set(data.event_codes) <= {1}
else {code: str(code) for code in data.event_codes}
)
if declared_causes != held_causes:
raise DataError(
f"the supplied LongitudinalData was built with causes={held_causes!r}, but "
f"this design declares causes={self.causes!r}; build the data from this "
"design, or correct the design"
)
[docs]
@dataclass(frozen=True)
class LongPanel:
"""Declare a longitudinal study read from a long panel: one row per unit and time.
Each row holds the covariates measured before that row's decision, the decision, and
whether the unit is still observed after it.
:meth:`cleverly.longitudinal.LongitudinalData.from_long_panel` states the layout, the
generated column names and the refusals. The estimands are those of
:class:`LongitudinalTreatment`, and one unit's rows are one observation.
Parameters
----------
unit : str
The unit identifier column.
time : str
The decision-time column.
treatment : str or sequence of str
The decision at each row. A list of columns is one vector node per time.
baseline : sequence of str
Unit-level covariates measured before the first decision. The list names at least one
column.
time_varying : sequence of str
Covariates measured at each row before that row's decision.
censoring : str or None
``1`` where the unit is still observed after the row's decision.
outcome : str or None
An end-of-study outcome, one value per unit.
event : str or mapping of str to str or None
A survival event indicator on each row, or one column per competing cause.
Declare exactly one of ``outcome`` and ``event``.
times : sequence or None
Every decision time, in time order. ``None`` uses the times the panel holds, sorted,
and needs numbers, dates or durations.
cluster : str or None
Independent-cluster identifier, constant within a unit.
weights : str or None
Unit-level probability-weight column.
weights_type : {"probability"}
Interpretation of ``weights``.
weights_estimated : bool
Whether the supplied weights were estimated from these data.
outcome_family : {"auto", "gaussian", "binomial"}
Outcome family used by the sequential regressions.
continuous_treatment : bool
Whether the treatment holds a continuous dose at every node.
treatment_support : Support or None
The declared conditional support of the dose at every node.
dose_atoms : sequence of float
The doses with a recorded point mass at every node.
treatment_probabilities : mapping of level to str, or None
The known treatment mechanism: each level mapped to the column that holds
``P(A_t = level | past)`` on each row. A tuple of ``(level, column)`` pairs, the form
the design stores, is also accepted.
censoring_probabilities : str or None
The column that holds ``P(C_t = 1 | past)`` on each row.
See Also
--------
LongitudinalTreatment : The same declaration on a wide frame.
CausalStudy : What a design is handed to along with the data.
cleverly.longitudinal.LongitudinalData.from_long_panel : The container it builds.
Examples
--------
>>> from cleverly import LongPanel
>>> design = LongPanel(
... unit="customer",
... time="week",
... treatment="offer",
... baseline=["age"],
... time_varying=["engagement"],
... censoring="active",
... outcome="bought",
... )
>>> design.time_varying
('engagement',)
"""
unit: str
time: str
treatment: str | Sequence[str]
baseline: Sequence[str]
time_varying: Sequence[str] = ()
censoring: str | None = None
outcome: str | None = None
event: str | Mapping[str, str] | None = None
times: Sequence[Any] | None = None
cluster: str | None = None
weights: str | None = None
weights_type: Literal["probability"] = "probability"
weights_estimated: bool = False
outcome_family: Family = "auto"
continuous_treatment: bool = False
treatment_probabilities: Mapping[Any, str] | None = None
censoring_probabilities: str | None = None
treatment_support: Any = None
dose_atoms: Sequence[float] = ()
def __post_init__(self) -> None:
object.__setattr__(self, "dose_atoms", tuple(float(atom) for atom in self.dose_atoms))
if (self.outcome is None) == (self.event is None):
raise DataError(
"LongPanel declares exactly one of outcome= (one end-of-study value per unit) "
"and event= (an event indicator on each row)"
)
if not isinstance(self.treatment, str):
object.__setattr__(self, "treatment", tuple(self.treatment))
object.__setattr__(self, "baseline", tuple(self.baseline))
object.__setattr__(self, "time_varying", tuple(self.time_varying))
if isinstance(self.event, Mapping):
object.__setattr__(self, "event", tuple(self.event.items()))
if self.times is not None:
object.__setattr__(self, "times", tuple(self.times))
if self.treatment_probabilities is not None:
declared = _stored_pairs(self.treatment_probabilities)
object.__setattr__(self, "treatment_probabilities", tuple(declared.items()))
@property
def outcome_role(self) -> Any:
"""The outcome the identified functional names: the outcome or the event columns."""
if self.outcome is not None:
return self.outcome
return dict(self.event) if isinstance(self.event, tuple) else self.event
[docs]
def prepare(self, data: Any) -> LongitudinalData:
"""Validate a long panel and construct its trajectory representation.
Parameters
----------
data : dataframe or LongitudinalData
A pandas or polars long panel, or data that this design prepared.
Returns
-------
LongitudinalData
Validated time-ordered representation.
"""
if isinstance(data, LongitudinalData):
self._check_prepared(data)
return data
built = LongitudinalData.from_long_panel(
data,
unit=self.unit,
time=self.time,
treatment=self.treatment,
baseline=self.baseline,
time_varying=self.time_varying,
censoring=self.censoring,
outcome=self.outcome,
event=dict(self.event) if isinstance(self.event, tuple) else self.event,
times=self.times,
cluster=self.cluster,
weights=self.weights,
weights_type=self.weights_type,
weights_estimated=self.weights_estimated,
family=self.outcome_family,
continuous_treatment=self.continuous_treatment,
treatment_probabilities=(
None if self.treatment_probabilities is None else dict(self.treatment_probabilities)
),
censoring_probabilities=self.censoring_probabilities,
)
return _declare_every_node(built, self.treatment_support, self.dose_atoms)
def _check_prepared(self, data: LongitudinalData) -> None:
"""Reconcile an already-built container with the roles this design declares."""
layout = data.panel
if layout is None:
raise DataError(
"the supplied LongitudinalData was built from a wide frame, but this design "
"declares a long panel; pass the panel, or use LongitudinalTreatment"
)
role: Any = None if self.event is None else self.outcome_role
event = () if role is None else (role,) if isinstance(role, str) else tuple(role.values())
expected: tuple[tuple[str, Any, Any], ...] = (
("unit", self.unit, layout.unit),
("time", self.time, layout.time),
(
"treatment",
(self.treatment,) if isinstance(self.treatment, str) else self.treatment,
layout.treatment,
),
("baseline", tuple(self.baseline), layout.baseline),
("time_varying", tuple(self.time_varying), layout.time_varying),
("censoring", self.censoring, layout.censoring),
("outcome", self.outcome, layout.outcome),
("event", event, layout.event),
("cluster", self.cluster, layout.cluster),
("weights", self.weights, layout.weights),
(
"treatment_probabilities",
None
if self.treatment_probabilities is None
else dict(self.treatment_probabilities),
None
if layout.treatment_probabilities is None
else dict(layout.treatment_probabilities),
),
(
"censoring_probabilities",
self.censoring_probabilities,
layout.censoring_probabilities,
),
("weights_type", self.weights_type, data.weight_spec.kind),
("weights_estimated", self.weights_estimated, data.weight_spec.estimated),
("continuous_treatment", self.continuous_treatment, data.is_continuous_node(1)),
("treatment_support", self.treatment_support, data.node_support(1)),
("dose_atoms", tuple(sorted(self.dose_atoms)), data.node_atoms(1)),
)
if self.times is not None:
expected += (("times", tuple(self.times), layout.decision_times),)
if self.outcome_family != "auto":
expected += (("outcome_family", self.outcome_family, data.family),)
for role, declared, held in expected:
if declared != held:
raise DataError(
f"the supplied LongitudinalData was built with {role}={held!r}, but this "
f"design declares {role}={declared!r}; build the data from this design, "
"or correct the design"
)
def _continuous_names(declared: Any) -> tuple[str, ...]:
"""The node names of a normalized ``continuous_treatment``: names, or (name, support)."""
return tuple(item[0] if isinstance(item, tuple) else item for item in declared)
def _named_declarations(declared: Any) -> dict[str, Any]:
"""The declared supports of a normalized ``continuous_treatment``, without ``None``."""
return {
item[0]: item[1] for item in declared if isinstance(item, tuple) and item[1] is not None
}
def _declared_continuous(declared: Any) -> Sequence[str] | Mapping[str, Any]:
"""A normalized ``continuous_treatment`` in the form ``from_frame`` reads."""
if declared and isinstance(declared[0], tuple):
return dict(declared)
return declared
def _point_atoms(declared: Any) -> Any:
"""Declared point atoms in one form: sorted floats, or sorted ``(component, atoms)`` pairs.
A vector treatment declares its atoms by component, as a mapping or as the pairs
:attr:`cleverly.data.CausalData.dose_atoms` holds.
"""
if isinstance(declared, Mapping):
declared = tuple(declared.items())
pairs = [entry for entry in declared if isinstance(entry, tuple)]
if pairs:
return tuple(
sorted((str(name), tuple(sorted(float(a) for a in atoms))) for name, atoms in pairs)
)
return tuple(sorted(float(atom) for atom in declared))
def _declared_atoms(declared: Any) -> Any:
"""The atoms as ``with_dose_declarations`` takes them: a mapping for a vector."""
if declared and isinstance(declared[0], tuple):
return dict(declared)
return declared
def _point_support(declared: Any) -> Any:
"""A point support for comparison: a vector's component supports become a dict."""
return dict(declared) if isinstance(declared, Mapping | tuple) else declared
def _held_supports(data: LongitudinalData) -> dict[str, Any]:
"""The declared support of each continuous node or component of ``data``, by name."""
if data.treatment_held:
support = data.node_support(1)
return {} if support is None else {data.treatment_names[0]: support}
held = {
name: data.node_support(time)
for time, name in enumerate(data.treatment_names, start=1)
if data.node_support(time) is not None
}
for time in range(1, data.n_times + 1):
held.update(data.component_supports(time))
return held
def _held_atoms(data: LongitudinalData) -> dict[str, tuple[float, ...]]:
"""The declared atoms of each continuous node or component of ``data``, by name."""
if data.treatment_held:
atoms = data.node_atoms(1)
return {data.treatment_names[0]: atoms} if atoms else {}
held = {
name: data.node_atoms(time)
for time, name in enumerate(data.treatment_names, start=1)
if data.node_atoms(time)
}
for time in range(1, data.n_times + 1):
held.update(data.component_atoms(time))
return held
def _declare_every_node(
data: LongitudinalData, support: Any, atoms: Sequence[float]
) -> LongitudinalData:
"""Declare one support and one atom set at every treatment node of ``data``."""
if support is None and not atoms:
return data
names = tuple(dict.fromkeys(data.treatment_names))
return data.with_dose_declarations(
{} if support is None else dict.fromkeys(names, support),
{} if not atoms else dict.fromkeys(names, atoms),
)
def _default_grid(data: LongitudinalData) -> tuple[float, ...]:
"""The grid ``grid=None`` gives on this container's times: 1 to the largest event time."""
times, codes = data.to_time_to_event()
events = times[codes > 0]
last = float(events.max()) if events.size else 0.0
if not last.is_integer():
return ()
return tuple(float(g) for g in range(1, int(last) + 1))
[docs]
@dataclass(frozen=True)
class ATE:
"""Compare population counterfactual means against one reference arm.
Parameters
----------
reference : Any or None
Treatment level used as the reference. ``None`` uses the design default.
See Also
--------
ATT : The same contrast restricted to the units in each comparison arm.
ATC : The same contrast restricted to the units in the reference arm.
RiskRatio : The same comparison on the ratio scale.
CounterfactualMean : The arm means the contrast is taken between.
Examples
--------
>>> from cleverly import ATE
>>> estimand = ATE(reference=0)
>>> estimand.reference
0
"""
reference: Any = None
name: str = field(default="ate", init=False)
@property
def definition(self) -> str:
"""Return the population-average contrast definition."""
return "average treatment effect, E[Y^a] - E[Y^reference]"
[docs]
@dataclass(frozen=True)
class ATT:
"""Compare counterfactual means among units in each comparison arm.
Parameters
----------
reference : Any or None
Treatment level used as the reference. ``None`` uses the design default.
See Also
--------
ATE : The same contrast over the whole population.
ATC : The same contrast restricted to the units in the reference arm.
cleverly.ControlledDirectEffect : An arm contrast holding an intermediate fixed.
Examples
--------
>>> from cleverly import ATT
>>> estimand = ATT(reference=0)
>>> estimand.reference
0
"""
reference: Any = None
name: str = field(default="att", init=False)
@property
def definition(self) -> str:
"""Return the treated-population contrast definition."""
return "average treatment effect among units receiving each comparison arm"
[docs]
@dataclass(frozen=True)
class ATC:
"""Compare counterfactual means among units in the reference arm.
Parameters
----------
reference : Any or None
Treatment level that defines the target population and comparison reference.
See Also
--------
ATE : The same contrast over the whole population.
ATT : The same contrast restricted to the units in each comparison arm.
CounterfactualMean : The arm means the contrast is taken between.
Examples
--------
>>> from cleverly import ATC
>>> estimand = ATC(reference=0)
>>> estimand.reference
0
"""
reference: Any = None
name: str = field(default="atc", init=False)
@property
def definition(self) -> str:
"""Return the reference-population contrast definition."""
return "average treatment effect among units receiving the reference arm"
[docs]
@dataclass(frozen=True)
class CounterfactualMean:
"""Request counterfactual outcome means by treatment level.
Parameters
----------
treatment : Any or None
One treatment level to retain. ``None`` reports every supported level.
See Also
--------
ATE : The contrast between these means.
RegimeMean : Means under declared regimens rather than fixed levels.
ModifiedTreatmentPolicy : Means under modified treatment policies.
Examples
--------
>>> from cleverly import CounterfactualMean
>>> estimand = CounterfactualMean(treatment=1)
>>> estimand.treatment
1
"""
treatment: Any = None
name: str = field(default="ey", init=False)
@property
def definition(self) -> str:
"""Return the counterfactual-mean definition."""
return (
"counterfactual mean under each treatment, E[Y^a]"
if self.treatment is None
else f"counterfactual mean under treatment {self.treatment!r}, E[Y^a]"
)
[docs]
@dataclass(frozen=True)
class NaturalCourseMean:
"""Request the observed-course outcome mean.
With complete outcomes this is the empirical mean of ``Y``. It reads the observed
law alone, so it estimates no nuisance and makes no exchangeability assumption.
When some outcome is actually missing and ``missingness=`` is declared, the same
target is ``E[m(A, W)]``, where ``m(A, W) = E(Y | Delta = 1, A, W)``. The
complete-case mean is a different parameter. That form estimates the outcome
regression and the response mechanism, and still estimates no treatment mechanism.
It needs missingness at random given ``(A, W)`` and response positivity. See
:doc:`observed-data extensions </technical-reference/point-treatment-tmle>` for the
supported compositions. The in-sample estimator accepts a binary outcome or a bounded
continuous outcome with declared ``q_bounds``. The cross-fitted estimator accepts a
binary outcome and one package-generated unstratified V-fold partition with at least
two folds. It uses one pooled fluctuation and evaluates the point and influence curve
on all stacked held-out rows. One fold would be the in-sample estimator, so the
cross-fitted estimator refuses it rather than report that estimate under the stacked
contract and its second-moment covariance rule. The in-sample estimator is ordinary
TMLE. The cross-fitted estimator is stacked CV-TMLE. Both are scalar estimators with
binary treatment. A fit outside these two estimators raises ``CapabilityError``
before any learner is fitted. A method declaration that is invalid on its own raises
``MethodConfigurationError`` when it is constructed. A sample without enough
respondents or nonrespondents for cross-fitting raises ``DataError`` before any
learner is fitted. The
:doc:`scope and refusals </technical-reference/scope-and-refusals>` page lists each
refusal, including joint targets, bootstrap inference, and linear or one-step
targeting.
Which of the two forms runs is settled by the data, not by the declaration. An
observation indicator that is one on every row leaves nothing missing, so the fit
takes the complete-outcome branch above and the narrow boundary does not apply to
it. The two agree there: with no missing outcome the empirical mean *is*
``E[m(A, W)]``. That branch is not the stacked estimator, so it follows the ordinary
cross-fitting rules rather than the contract above. Those rules include the outcome
scale: a cross-fitted fit of a continuous outcome needs a declared ``q_bounds``, or
``cross_fit=False``.
"""
name: str = field(default="ey_obs", init=False)
definition: str = field(default="natural-course outcome mean, E[Y]", init=False)
[docs]
@dataclass(frozen=True)
class PopulationAttributableRisk:
"""Request the natural-course risk minus a reference counterfactual risk.
Parameters
----------
reference : Any or None
Treatment level representing removal of the exposure.
"""
reference: Any = None
name: str = field(default="par", init=False)
definition: str = field(default="E[Y] - E[Y^reference]", init=False)
[docs]
@dataclass(frozen=True)
class PopulationAttributableFraction:
"""Request the preventable fraction relative to the natural-course risk.
Parameters
----------
reference : Any or None
Treatment level representing removal of the exposure.
"""
reference: Any = None
name: str = field(default="paf", init=False)
definition: str = field(default="1 - E[Y^reference] / E[Y]", init=False)
[docs]
@dataclass(frozen=True)
class RiskRatio:
"""Compare counterfactual risks on the ratio scale.
Parameters
----------
reference : Any or None
Denominator treatment level. ``None`` uses the design default.
See Also
--------
OddsRatio : The same comparison on the odds scale.
ATE : The same comparison on the difference scale.
cleverly.sensitivity.evalue.evalue : How much confounding would explain a ratio away.
Examples
--------
>>> from cleverly import RiskRatio
>>> estimand = RiskRatio(reference=0)
>>> estimand.reference
0
"""
reference: Any = None
name: str = field(default="rr", init=False)
definition: str = field(default="counterfactual risk ratio", init=False)
[docs]
@dataclass(frozen=True)
class OddsRatio:
"""Compare counterfactual odds for a binary outcome.
Parameters
----------
reference : Any or None
Denominator treatment level. ``None`` uses the design default.
"""
reference: Any = None
name: str = field(default="or", init=False)
definition: str = field(default="counterfactual odds ratio", init=False)
[docs]
@dataclass(frozen=True)
class RegimeMean:
"""Request mean outcomes under declared treatment regimens.
Parameters
----------
regimens : Any
Point-treatment interventions or longitudinal regimens. On a point-treatment
design, a sequence of treatment levels and regimes, such as ``(0, Static(1))``.
``CausalStudy.identify`` refuses a ``Shift`` or an ``Incremental`` here.
reference : str or None
Regimen label used as the reference for related contrasts.
horizons : sequence of int or None
Follow-up times to report for a longitudinal outcome.
"""
regimens: Any
reference: str | None = None
horizons: Sequence[int] | None = None
name: str = field(default="ey_regime", init=False)
definition: str = field(default="mean outcome under each declared regime", init=False)
[docs]
@dataclass(frozen=True)
class RegimeContrast:
"""Compare each declared regimen against one reference regimen.
Parameters
----------
regimens : Any
Point-treatment interventions or longitudinal regimens. On a point-treatment
design, a sequence of treatment levels and regimes, such as ``(0, Static(1))``.
``CausalStudy.identify`` refuses a ``Shift`` or an ``Incremental`` here.
reference : str or None
Reference regimen label. ``None`` uses the first regimen.
horizons : sequence of int or None
Follow-up times to report for a longitudinal outcome.
See Also
--------
RegimeMean : The regimen means this contrast is taken between.
ATE : The same comparison when treatment is given once.
cleverly.LongitudinalTreatment : The design a regimen is declared against.
Examples
--------
Two static longitudinal regimens, compared against never treating:
>>> from cleverly import RegimeContrast
>>> estimand = RegimeContrast(
... regimens={"always": 1, "never": 0}, reference="never"
... )
>>> estimand.reference
'never'
"""
regimens: Any
reference: str | None = None
horizons: Sequence[int] | None = None
name: str = field(default="ate_regime", init=False)
definition: str = field(default="contrast of each regime against the reference", init=False)
[docs]
@dataclass(frozen=True)
class LearnedRuleValue:
"""Request the fold-average value of a treatment rule learned inside each training fold.
The target is data-adaptive. It is the average over the outer folds of the value of
the rule learned on each fold's training rows, so it depends on the realized split and
on the fitted rules. It is not the value of one rule fitted on all rows, and it is not
the value of the optimal rule. Van der Laan and Luedtke (2015), Section 7 and
Appendix B, define it and its cross-validated TMLE. The interval can under-cover when
the effect is zero or near zero for a share of units, where no limiting rule exists.
The fit needs the fold-evaluated CV-TMLE, ``CrossFitting(enabled=True,
fold_evaluation=True)``. The default ``CrossFitting()`` refuses before any learner
and names that remedy.
Parameters
----------
rule : LearnedRule
The rule class, and the label the reported parameter carries.
See Also
--------
cleverly.interventions.LearnedRule : The rule this estimand holds.
RegimeMean : The value of a known rule, fixed before the fit.
cleverly.CrossFitting : The fold-evaluated configuration this estimand needs.
Examples
--------
>>> from cleverly import LearnedRuleValue
>>> estimand = LearnedRuleValue()
>>> estimand.name, estimand.rule.name
('ey_learned_rule', 'learned rule')
"""
rule: LearnedRule = field(default_factory=LearnedRule)
name: str = field(default="ey_learned_rule", init=False)
definition: str = field(
default=(
"average over the outer folds of the value of the rule learned on each fold's "
"training rows"
),
init=False,
)
def __post_init__(self) -> None:
if not isinstance(self.rule, LearnedRule):
raise DataError(
"LearnedRuleValue(rule=...) takes a cleverly.interventions.LearnedRule, not "
f"{self.rule!r}"
)
[docs]
@dataclass(frozen=True)
class ModifiedTreatmentPolicy:
"""Request mean outcomes under modified treatment policies.
Parameters
----------
policies : sequence of Policy
Named policies to evaluate: ``Shift``, ``Scale``, ``Piecewise`` or
``ModifiedPolicy``. A ``RiskRatioTilt`` reports its own estimand, ``ey_rr_tilt``, which
has no typed request; this class refuses it.
reference : str or None
Policy label retained as the comparison reference.
See Also
--------
ModifiedTreatmentPolicyEffect : The contrast between these means.
cleverly.interventions.Shift : The additive dose shift, the simplest policy.
cleverly.interventions.check_policy_support : Whether the shifted dose is supported.
Examples
--------
Raise every dose by one unit without imposing a cap:
>>> from cleverly import ModifiedTreatmentPolicy
>>> from cleverly.interventions import Shift
>>> estimand = ModifiedTreatmentPolicy(policies=[Shift(delta=1.0, cap=None, name="up1")])
>>> estimand.policies[0].cap is None
True
"""
policies: Sequence[Policy]
reference: str | None = None
name: str = field(default="ey_policy", init=False)
definition: str = field(default="mean outcome under each modified treatment policy", init=False)
def __post_init__(self) -> None:
_refuse_typed_tilts(self.policies, type(self).__name__)
[docs]
@dataclass(frozen=True)
class ModifiedTreatmentPolicyEffect:
"""Compare modified treatment policies.
Parameters
----------
policies : sequence of Policy
Named policies to compare. A ``RiskRatioTilt`` is refused, as in
:class:`ModifiedTreatmentPolicy`.
reference : str or None
Reference policy label. ``None`` uses the first policy.
"""
policies: Sequence[Policy]
reference: str | None = None
name: str = field(default="ate_policy", init=False)
definition: str = field(default="contrast of modified treatment policies", init=False)
def __post_init__(self) -> None:
_refuse_typed_tilts(self.policies, type(self).__name__)
def _refuse_typed_tilts(policies: Sequence[Any], holder: str) -> None:
"""Refuse a risk-ratio tilt in a typed policy request, which reports ``ey_policy``."""
from .interventions import RiskRatioTilt
tilts = [getattr(item, "name", "") for item in policies if isinstance(item, RiskRatioTilt)]
if tilts:
raise CapabilityError(
f"{holder} holds the risk-ratio tilt(s) {tilts}. A tilt reports its own estimands, "
"ey_rr_tilt and ate_rr_tilt, and no typed study request reports them.",
reason=RefusalReason.API_REDIRECTION,
alternative=(
"Fit TMLE(policies=[RiskRatioTilt(1.0), RiskRatioTilt(delta)]) directly and "
"read estimates['ey_rr_tilt[...]']."
),
)
[docs]
@dataclass(frozen=True)
class IncrementalMean:
"""Request means under incremental propensity-score interventions.
Parameters
----------
interventions : sequence of Incremental
Odds multipliers to evaluate.
reference : str or None
Intervention label retained as the comparison reference.
horizons : sequence of int or None
Follow-up times to report for a longitudinal outcome. On a repeated outcome, each
horizon is the end-of-study fit of the panel truncated after it. A
point-treatment design refuses it.
"""
interventions: Sequence[Incremental]
reference: str | None = None
horizons: Sequence[int] | None = None
name: str = field(default="ey_ipsi", init=False)
definition: str = field(
default="mean under each incremental propensity intervention", init=False
)
[docs]
@dataclass(frozen=True)
class IncrementalEffect:
"""Compare incremental propensity-score interventions.
Parameters
----------
interventions : sequence of Incremental
Odds multipliers to compare.
reference : str or None
Reference intervention label. ``None`` uses the first intervention.
horizons : sequence of int or None
Follow-up times to report for a longitudinal outcome, as for
:class:`IncrementalMean`.
See Also
--------
IncrementalMean : The intervention means this contrast is taken between.
cleverly.interventions.Incremental : The odds multiplier an intervention applies.
ATE : The same comparison for a treatment set to a fixed level.
Examples
--------
Double the odds of treatment, and compare that against halving them:
>>> from cleverly import IncrementalEffect
>>> from cleverly.interventions import Incremental
>>> estimand = IncrementalEffect(
... interventions=[
... Incremental(delta=0.5, name="halved"),
... Incremental(delta=2.0, name="doubled"),
... ],
... reference="halved",
... )
>>> estimand.reference
'halved'
"""
interventions: Sequence[Incremental]
reference: str | None = None
horizons: Sequence[int] | None = None
name: str = field(default="ate_ipsi", init=False)
definition: str = field(default="contrast of incremental propensity interventions", init=False)
[docs]
@dataclass(frozen=True)
class MSMProjection:
"""Project counterfactual means onto a marginal structural model.
Parameters
----------
model : MSM
Working model that defines the projection terms and weights.
regimens : Any or None
Longitudinal regimen cells. Leave unset for a point-treatment projection.
horizons : sequence of int or None
Follow-up times included in a longitudinal projection.
"""
model: MSM
regimens: Any = None
horizons: Sequence[int] | None = None
name: str = field(default="msm", init=False)
definition: str = field(
default="projection of counterfactual means onto a working model", init=False
)
[docs]
@dataclass(frozen=True)
class ControlledDirectEffect:
"""Compare treatments while fixing a binary intermediate variable.
Parameters
----------
intermediate : float
Fixed intermediate level. The evidenced implementation accepts zero or one.
contrast : ATE, ATT, ATC, RiskRatio, or OddsRatio
Arm contrast evaluated at the fixed intermediate level.
See Also
--------
ATE : The same contrast with nothing held fixed.
cleverly.PointTreatment : Where the ``intermediate`` column is declared.
cleverly.datasets.make_cde : A process with a known controlled direct effect.
Examples
--------
>>> from cleverly import ATE, ControlledDirectEffect
>>> estimand = ControlledDirectEffect(intermediate=0.0, contrast=ATE())
>>> estimand.intermediate
0.0
"""
intermediate: float
contrast: ATE | ATT | ATC | RiskRatio | OddsRatio = field(default_factory=ATE)
name: str = field(default="controlled_direct_effect", init=False)
def __post_init__(self) -> None:
# Identification reads ``self.contrast.name`` before anything else, so a string here
# would surface as ``AttributeError: 'str' object has no attribute 'name'`` from deep
# inside the provider rather than as a refusal naming what to pass instead.
if not isinstance(self.contrast, (ATE, ATT, ATC, RiskRatio, OddsRatio)):
raise DataError(
"ControlledDirectEffect(contrast=...) takes a typed arm contrast -- "
f"ATE(), ATT(), ATC(), RiskRatio(), or OddsRatio() -- not {self.contrast!r}"
)
@property
def definition(self) -> str:
"""Return the controlled-direct-effect definition."""
return f"controlled direct effect with the intermediate fixed at {self.intermediate}"
PointEstimand = (
ATE
| ATT
| ATC
| CounterfactualMean
| NaturalCourseMean
| PopulationAttributableRisk
| PopulationAttributableFraction
| RiskRatio
| OddsRatio
| RegimeMean
| RegimeContrast
| LearnedRuleValue
| ModifiedTreatmentPolicy
| ModifiedTreatmentPolicyEffect
| IncrementalMean
| IncrementalEffect
| MSMProjection
| ControlledDirectEffect
)
#: The typed object to reach for when a caller passes a string spelling. Only the
#: names ``TMLE(estimands=(...,))`` accepts as strings are listed; anything else gets the
#: generic pointer to the roster in ``cleverly.__all__``.
_STRING_ESTIMANDS: dict[str, str] = {
"ate": "ATE()",
"att": "ATT()",
"atc": "ATC()",
"ey": "CounterfactualMean()",
"ey1": "CounterfactualMean(treatment=1)",
"ey0": "CounterfactualMean(treatment=0)",
"ey_obs": "NaturalCourseMean()",
"par": "PopulationAttributableRisk()",
"paf": "PopulationAttributableFraction()",
"rr": "RiskRatio()",
"or": "OddsRatio()",
"ey_regime": "RegimeMean(regimens=...)",
"ate_regime": "RegimeContrast(regimens=...)",
"ey_learned_rule": "LearnedRuleValue()",
"ey_policy": "ModifiedTreatmentPolicy(policies=...)",
"ate_policy": "ModifiedTreatmentPolicyEffect(policies=...)",
"ey_ipsi": "IncrementalMean(interventions=...)",
"ate_ipsi": "IncrementalEffect(interventions=...)",
"msm": "MSMProjection(model=...)",
}
#: Why :class:`~cleverly.longitudinal.IPSI` is unavailable for a plan with a policy node.
IPSI_POLICY_REFUSAL = Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"IPSI estimates no density ratio for a modified treatment policy, a joint policy or a "
"risk-ratio tilt, so it has no estimator for a plan with such a node.",
alternative="Fit the plans by LTMLE or SDR.",
roadmap="F33",
)
def _narrow_bands(
result: Any, retained: Mapping[str, Any], method: EstimationMethod
) -> SimultaneousBands | None:
"""Joint bands over the parameters this effect reports, not the ones the engine computed.
A typed estimand can ask for a subset of what its engine produces: ``CounterfactualMean``
naming an arm of a multi-arm fit, or a longitudinal target whose engine reports means and
contrasts together. Filtering ``estimates`` alone left ``simultaneous`` behind, so the
result carried a critical value quantifying a family it no longer contained and bands
keyed to parameters it could not index -- and ``summary()`` printed them.
The bands are recomputed rather than dropped, from the same influence curves and the same
``alpha``, draw count, multiplier distribution, seed and cluster structure the engine used.
That makes them the bands the engine itself would have produced had it been asked for this
family alone. Below two parameters there is no joint question left, and both engines
already return ``None`` there rather than a one-parameter band.
A parameter that supplies no inference on its own, such as a constant-node risk of a
longitudinal fit, stays out of the family, as the engine's band leaves it out. The band is
a joint statement over the retained parameters that supply inference.
"""
bands = result.simultaneous
if bands is None:
return None
family = {name: estimate for name, estimate in retained.items() if estimate.supplies_inference}
if set(family) == set(bands.bands):
return bands
if len(family) < 2:
return None
narrowed = simultaneous_bands(
family,
alpha=bands.alpha,
n_replicates=bands.n_replicates,
kind=bands.kind,
random_state=getattr(method, "runtime", Runtime()).random_state,
cluster=result.data.cluster,
)
return replace(narrowed, support=combine_support([bands.support, narrowed.support]))
def _narrow_bootstrap(bootstrap: Any, retained: Mapping[str, Any]) -> Any:
"""Bootstrap draws for the parameters this effect reports, and no others.
Unlike a joint band, per-parameter draws stay correct under narrowing -- each estimand's
resampling distribution is its own. What they must not do is outlive their parameter: a
result narrowed to one arm would still hand out ``bootstrap.draws`` for the arms it
dropped, keyed by names it cannot itself index.
"""
if bootstrap is None:
return None
return replace(
bootstrap,
draws={name: values for name, values in bootstrap.draws.items() if name in retained},
)
[docs]
@dataclass(frozen=True)
class BackdoorMeanContrast:
"""Store a normalized observed-data functional for an estimator adapter.
Users normally receive this object through :meth:`CausalStudy.identify`.
Parameters
----------
outcome : Any
Outcome column declaration.
treatment : Any
Treatment column declaration.
adjustment : tuple of str
Adjustment columns or baseline history.
target : str
Registered engine target.
axis : str
Parameter axis, such as ``"arm"``, ``"regimen"``, or ``"policy"``.
reference : Any or None
Reference level or label.
interventions : Any
Interventions, shifts, or regimens evaluated by the target.
horizons : tuple of int or None
Requested longitudinal follow-up times.
msm : MSM or None
Marginal structural model for a projection target.
intermediate : float or None
Fixed intermediate level for a controlled direct effect.
longitudinal : bool
Whether the functional uses sequential identification.
missingness : str or None
Outcome-observation indicator used by the observed-data regression.
treatment_missingness : str or None
Treatment-observation indicator used by the observed-data regression.
intermediate_name : str or None
Intermediate column fixed by a controlled direct effect.
treatment_levels : tuple
Declared treatment levels, in the data container's stable order.
treatment_value : Any or None
One retained treatment level for a targeted counterfactual mean.
history : tuple of tuple
For a longitudinal functional, each treatment node with the set that its treatment
mechanism conditions on: the covariates measured before the node, then the earlier
treatment nodes. Empty for a point treatment.
"""
outcome: Any
treatment: Any
adjustment: tuple[str, ...]
target: str
axis: str = "arm"
reference: Any = None
interventions: Any = ()
horizons: tuple[int, ...] | None = None
msm: MSM | None = None
intermediate: float | None = None
longitudinal: bool = False
missingness: str | None = None
intermediate_name: str | None = None
treatment_levels: tuple[Any, ...] = ()
treatment_missingness: str | None = None
treatment_value: Any = None
history: tuple[tuple[str, tuple[str, ...]], ...] = ()
@property
def reference_arm(self) -> Any:
"""Return the arm this functional's contrasts are taken against.
One implementation, because the printed expression and the printed positivity
statement both resolve it and a reader who saw them disagree could not tell which
arm the fit used. ``"reference"`` stands in when the record declares no reference
and no :attr:`treatment_levels`, as on a continuous treatment, which leaves
nothing to name.
"""
if self.reference is not None:
return self.reference
if self.treatment_levels:
return self.treatment_levels[0]
return "reference"
@property
def selected_arm(self) -> Any:
"""Return the one arm an evaluated-mean functional reports, or ``None`` for all.
``ey1`` and ``ey0`` name an arm through the target rather than through
``treatment_value``, which is why the rule is written here rather than at each of
the three places that used to repeat it.
"""
if self.treatment_value is not None:
return self.treatment_value
if self.target in {"ey1", "ey0"} and self.treatment_levels:
return self.treatment_levels[1 if self.target == "ey1" else 0]
return None
@property
def expression(self) -> str:
"""Return a readable expression for the identified functional.
Raises
------
CapabilityError
When an arm-indexed target has no entry below. The alternative is to print
some other estimand's functional into the record that ``summary()`` shows and
``save()`` persists, which is a wrong answer rather than a missing one.
"""
if self.longitudinal:
return "sequential g-formula under the declared treatment regimen"
def regression(arm: Any = "a") -> str:
conditions = [f"{self.treatment}={arm}"]
if self.intermediate is not None:
name = self.intermediate_name or "Z"
conditions.append(f"{name}={self.intermediate:g}")
if self.treatment_missingness is not None:
conditions.append(f"{self.treatment_missingness}=1")
if self.missingness is not None:
conditions.append(f"{self.missingness}=1")
conditions.append("W")
return f"E({self.outcome} | {', '.join(conditions)})"
support = f" for a in {list(self.treatment_levels)!r}" if self.treatment_levels else ""
if self.axis == "learned_rule":
treated, control = self.treatment_levels[1], self.treatment_levels[0]
return (
f"(1/V) sum_v E_W[{regression('d_v(W)')}], where d_v(W) = 1 when "
f"{regression(repr(treated))} - {regression(repr(control))} > 0 as learned "
"on the training rows of outer fold v"
)
if self.axis == "arm":
reference = self.reference_arm
comparisons = tuple(level for level in self.treatment_levels if level != reference)
comparison_domain = f" for a in {list(comparisons)!r}" if comparisons else ""
def marginal_mean(arm: Any) -> str:
return f"E_W[{regression(arm)}]"
if self.target == "ey_obs":
if self.missingness is not None:
return (
f"E_{{{self.treatment},W}}[E({self.outcome} | "
f"{self.missingness}=1, {self.treatment}, W)]"
)
return f"E({self.outcome}) under the observed treatment course"
if self.target == "att":
return (
f"E_{{W | {self.treatment}=a}}[{regression('a')} - "
f"{regression(repr(reference))}]{comparison_domain}"
)
if self.target == "atc":
return (
f"E_{{W | {self.treatment}={reference!r}}}[{regression('a')} - "
f"{regression(repr(reference))}]{comparison_domain}"
)
if self.target in {"par", "paf"}:
counterfactual = marginal_mean(repr(reference))
observed = (
f"E_{{{self.treatment},W}}[E({self.outcome} | "
f"{self.missingness}=1, {self.treatment}, W)]"
if self.missingness is not None
else f"E({self.outcome})"
)
if self.target == "par":
return f"{observed} - {counterfactual}"
return f"1 - {counterfactual} / {observed}"
if self.target in {"ey", "ey1", "ey0"}:
selected = self.selected_arm
if selected is not None:
return marginal_mean(repr(selected))
return f"{marginal_mean('a')}{support}"
if self.target == "ate":
return f"{marginal_mean('a')} - {marginal_mean(repr(reference))}{comparison_domain}"
if self.target == "rr":
return f"{marginal_mean('a')} / {marginal_mean(repr(reference))}{comparison_domain}"
if self.target == "or":
mean_a = marginal_mean("a")
mean_r = marginal_mean(repr(reference))
return (
f"({mean_a} / (1 - {mean_a})) / ({mean_r} / (1 - {mean_r})){comparison_domain}"
)
raise CapabilityError(
f"the arm-indexed target {self.target!r} has no identified expression; add "
"one beside the estimands above rather than reporting the counterfactual "
"mean's functional under another estimand's name. This is a package defect; "
"report it.",
reason=RefusalReason.INVALID_DECLARATION,
)
return (
f"identified {self.axis}-indexed plug-in functional of {regression()}{support} "
"with influence correction"
)
[docs]
@dataclass(frozen=True)
class ParameterKey:
"""Describe the structured identity behind a parameter alias.
Parameters
----------
alias : str
Stable key used in result mappings.
estimand : str
Registered estimand name.
axis : str
Dimension indexed by ``value``.
value : Any or None
Arm, intervention, or cell value.
reference : Any or None
Reference value used by a contrast.
stratum : tuple or None
Target stratum values.
regimen : str or None
Longitudinal regimen label.
horizon : int or None
Follow-up time.
cause : str or None
Competing-risk cause label.
term : str or None
Marginal structural model term.
"""
alias: str
estimand: str
axis: str = "arm"
value: Any = None
reference: Any = None
stratum: tuple[Any, ...] | None = None
regimen: str | None = None
horizon: int | None = None
cause: str | None = None
term: str | None = None
def _point_target(estimand: PointEstimand, data: CausalData) -> tuple[Any, str]:
"""Return the engine estimand and target selected by one typed point question."""
actual = estimand.contrast if isinstance(estimand, ControlledDirectEffect) else estimand
target = actual.name
if (
isinstance(actual, CounterfactualMean)
and actual.treatment is not None
and data.is_binary_treatment
and actual.treatment in data.treatment_levels
):
code = float(data.treatment_levels.index(actual.treatment))
target = "ey1" if code == 1.0 else "ey0"
return actual, target
@dataclass(frozen=True)
class _ArmScope:
"""The arms one point-treatment functional evaluates, resolved once.
``_point_positivity`` and ``_point_mechanism_scope`` were the same ladder twice --
same unwrap of a controlled direct effect, same reference-arm default, same
selected-arm rule, same order of cases -- and they had already drifted apart. They
now share this resolution and differ only in the sentence each case writes.
Parameters
----------
treatment : str
The treatment column, as the design names it.
reference : Any
The arm contrasts are taken against.
selected : Any or None
The one arm an evaluated-mean functional reports, or ``None`` for all arms.
levels : tuple
The supported treatment levels, empty for a continuous dose.
continuous : bool
Whether the design declared a continuous dose rather than arms.
"""
treatment: str
reference: Any
selected: Any
levels: tuple[Any, ...]
continuous: bool
def _all_arm_statements(scope: _ArmScope) -> tuple[str, str]:
"""State support at every declared arm: the shared case for ``ate``, ``rr`` and ``or``."""
return (
f"positivity: P({scope.treatment} = a | W) > 0 almost surely for every "
f"supported treatment level a in {list(scope.levels)!r}",
f"for every supported treatment level b in {list(scope.levels)!r} "
"almost surely over the target population",
)
def _evaluated_arm_statements(scope: _ArmScope) -> tuple[str, str]:
"""State support at the one arm a counterfactual mean evaluates."""
if scope.selected is None:
return _all_arm_statements(scope)
return (
f"positivity for the evaluated arm: P({scope.treatment} = {scope.selected!r} | W) "
"> 0 almost surely",
f"at b = {scope.selected!r} almost surely over the target population",
)
def _treated_population_statements(scope: _ArmScope) -> tuple[str, str]:
"""State the reference-arm support the effect on the treated needs."""
return (
f"positivity for the reference arm within each comparison population: "
f"P({scope.treatment} = {scope.reference!r} | W) > 0 wherever "
f"P({scope.treatment} = a | W) > 0, for every comparison arm a",
f"for b in {{a, {scope.reference!r}}} wherever P({scope.treatment} = a | W) > 0",
)
def _reference_population_statements(scope: _ArmScope) -> tuple[str, str]:
"""State the comparison-arm support the effect on the controls needs."""
return (
f"positivity for each comparison arm within the reference population: "
f"P({scope.treatment} = a | W) > 0 wherever "
f"P({scope.treatment} = {scope.reference!r} | W) > 0, for every comparison arm a",
f"for b in {{a, {scope.reference!r}}} wherever "
f"P({scope.treatment} = {scope.reference!r} | W) > 0",
)
def _reference_intervention_statements(scope: _ArmScope) -> tuple[str, str]:
"""State the support an attributable effect needs at its reference intervention."""
return (
f"positivity for the reference intervention: P({scope.treatment} = "
f"{scope.reference!r} | W) > 0 almost surely",
_all_arm_statements(scope)[1],
)
#: One row per target that can reach :func:`_arm_scope`: the treatment support that
#: target's functional needs, and the arm and population cells a fitted intermediate or
#: response mechanism must cover. Every estimand's prose is its own auditable entry. A
#: target reaching :func:`_arm_scope` without a row is refused rather than given the
#: counterfactual mean's sentence under its own name.
#:
#: ``ey_obs`` is the one registered arm target with no row, and it needs none: its record
#: keeps only the complete-outcome assumption, so the two-arm positivity sentence that
#: resolves the arms is dropped before the branch that would resolve them. Its functional
#: asks nothing of the treatment mechanism, so it has no treatment support to state. A
#: row that returned an empty sentence instead was dropped in silence by the caller, which
#: is the shape :func:`_arm_scope` now refuses.
#: ``test_the_arm_statement_rows_cover_every_target_that_can_reach_them`` compares these
#: keys with the registry and pins that exclusion.
_ARM_STATEMENTS: dict[str, Callable[[_ArmScope], tuple[str, str]]] = {
"ate": _all_arm_statements,
"rr": _all_arm_statements,
"or": _all_arm_statements,
"msm": _all_arm_statements,
"ey": _evaluated_arm_statements,
"ey1": _evaluated_arm_statements,
"ey0": _evaluated_arm_statements,
"att": _treated_population_statements,
"atc": _reference_population_statements,
"par": _reference_intervention_statements,
"paf": _reference_intervention_statements,
}
_CONTINUOUS_DOSE_STATEMENTS = (
"positivity for the evaluated doses: the conditional treatment density is "
"positive wherever the target functional gives a dose positive weight",
"wherever the target functional gives a dose positive weight",
)
def _design_has_arms(design: PointTreatment) -> bool:
"""Whether the design's treatment has arms, read from :attr:`PointTreatment.support`.
Every axis reader asks this one function, and it names each support it knows. A
support it does not know raises, so a new kind of treatment cannot take the arm branch
of a reader written before it existed.
"""
support = design.support
if support == "discrete":
return True
if support in ("continuous", "vector"):
return False
raise CleverlyError(
f"treatment support {support!r} has no reading in the study spine; the supports "
"are 'discrete', 'continuous' and 'vector'"
)
def _vector_statements(design: PointTreatment) -> tuple[str, str]:
"""The positivity sentence and the mechanism scope of a vector with a dose."""
columns = ", ".join(design.treatment_columns)
return (
f"positivity for the evaluated treatment values: the joint conditional law of "
f"({columns}) given W has a positive density wherever the target functional gives "
"a value positive weight",
"wherever the target functional gives a value positive weight",
)
def _arm_scope(
design: PointTreatment, functional: BackdoorMeanContrast
) -> tuple[_ArmScope, str, str]:
"""Resolve the arms once, and return the two statements the target's row writes.
Returns
-------
tuple
The resolved scope, the positivity sentence, and the mechanism-scope phrase.
The positivity sentence is empty for a target that needs none.
Raises
------
CapabilityError
When the target has no row in :data:`_ARM_STATEMENTS`, or when its row states no
treatment support. An empty sentence is not a statement that the target needs
none: the caller cannot tell it from a row nobody wrote, and appending it would
put a blank assumption in the record.
"""
scope = _ArmScope(
treatment=design.treatment_name,
reference=functional.reference_arm,
selected=functional.selected_arm,
levels=functional.treatment_levels,
continuous=not _design_has_arms(design),
)
if design.support == "vector":
return (scope, *_vector_statements(design))
if scope.continuous:
return (scope, *_CONTINUOUS_DOSE_STATEMENTS)
row = _ARM_STATEMENTS.get(functional.target)
if row is None:
raise CapabilityError(
f"the target {functional.target!r} declares the two-arm positivity assumption "
"but states no design-bound support of its own; add its row beside the "
"estimands in _ARM_STATEMENTS rather than printing another estimand's "
"positivity statement under its name. This is a package defect; report it.",
reason=RefusalReason.INVALID_DECLARATION,
)
positivity, cells = row(scope)
if not positivity:
raise CapabilityError(
f"the row for {functional.target!r} in _ARM_STATEMENTS states no treatment "
"support; a target that reaches here declares the two-arm positivity "
"assumption, so write the support its functional needs rather than an empty "
"sentence the record would drop without saying so. This is a package defect; "
"report it.",
reason=RefusalReason.INVALID_DECLARATION,
)
return (scope, positivity, cells)
def _conditioning_event_statement(design: PointTreatment, scope: _ArmScope, target: str) -> str:
"""State which conditioning event a conditional effect needs to be defined.
Raises
------
CapabilityError
For a target that declares the conditioning-event assumption but names no
population of its own.
"""
if target == "att":
return (
f"the conditioning event has positive probability: "
f"P({design.treatment_name} = a) > 0 for every comparison arm a"
)
if target == "atc":
return (
f"the reference conditioning event has positive probability: "
f"P({design.treatment_name} = {scope.reference!r}) > 0"
)
raise CapabilityError(
f"the target {target!r} declares a conditioning event but names no population "
"it conditions on; state which arm's units carry the parameter rather than "
"keeping another estimand's sentence. This is a package defect; report it.",
reason=RefusalReason.INVALID_DECLARATION,
)
def _bound_dr_condition(
base: Identification,
design: PointTreatment,
functional: BackdoorMeanContrast,
*,
direct: bool,
) -> str:
"""State the double-robustness condition this design leaves, in one sentence.
A design decides which mechanisms enter the remainder, and the registry cannot: a
declared response mechanism or a fixed intermediate makes the mechanism half a
*product*, so reporting "either Qbar or g" would name a fit where g alone is
sufficient when it is not.
Only the opening clause is the design's to write. A registry that states the
mechanism half as a condition carries
:data:`~cleverly.targets.builtin.DR_DESIGN_CONDITIONAL_CLAUSE`, and the design-bound
statement resolves it, so the whole sentence goes. Anything else a registered
sentence says is a fact about the estimand: the extra term ``att`` and ``atc`` carry
for the randomness of the conditioning event, and the note ``msm`` carries that beta
is a projection whatever the working model does. Those tails follow
:data:`~cleverly.targets.builtin.DR_MECHANISM_HEAD`, so the head is replaced and the
tail is kept. Replacing the whole string deleted the target-specific half of every
conditional and projection record under a ``missingness=`` design. A registered
sentence carrying neither is extended rather than rewritten, because a study that
cannot find the head it would replace still must not delete the sentence.
Parameters
----------
base : Identification
The registered theory for the target.
design : PointTreatment
The declared column roles and mechanisms.
functional : BackdoorMeanContrast
The design-bound functional the study identified.
direct : bool
Whether the design fixes an intermediate level.
Returns
-------
str
The condition under which this fit's estimate stays consistent.
"""
target = functional.target
missingness = functional.missingness
if target == "ey_obs" and missingness is None:
return "the complete-data empirical mean has no nuisance-model remainder"
if target == "ey_obs":
return (
"consistent if either m(A, W) = E(Y | Delta = 1, A, W) is consistent "
f"or P({missingness} = 1 | {design.treatment_name}, W) is consistent; no "
"treatment mechanism enters the remainder"
)
if target in POPULATION_INTERVENTION_TARGETS and missingness is not None:
treatment, reference = design.treatment_name, functional.reference_arm
return (
f"consistent if m({treatment}, W) = E({design.outcome} | {missingness} = 1, "
f"{treatment}, W) is consistent, or if both P({missingness} = 1 | {treatment}, W) "
f"and the product P({treatment} = {reference!r} | W) P({missingness} = 1 | "
f"{treatment} = {reference!r}, W) are consistent; a correct product with a wrong "
"response mechanism does not rescue the natural-course term"
)
if direct:
mechanism = "g * q_z" + (" * pi" if missingness is not None else "")
bound = (
"consistent if either Qbar(A, z, W) is consistent or the full mechanism "
f"product {mechanism} is consistent"
)
elif (
missingness is not None
and functional.axis in {"arm", "msm"}
and target not in POPULATION_INTERVENTION_TARGETS
):
bound = (
"consistent if either Qbar(A, W) is consistent or the full mechanism product "
f"g * P({missingness} = 1 | {design.treatment_name}, W) is consistent; a consistent "
"treatment mechanism alone is not sufficient once outcomes are missing"
)
else:
return base.dr_condition
if base.dr_condition.endswith(DR_DESIGN_CONDITIONAL_CLAUSE):
return bound
kept = base.dr_condition.removeprefix(DR_MECHANISM_HEAD)
if kept != base.dr_condition:
return bound + kept
return f"{base.dr_condition}. Under this design, {bound}"
def _point_identification(
base: Identification,
design: PointTreatment,
functional: BackdoorMeanContrast,
) -> Identification:
"""Bind registered theory to the treatment levels and active design mechanisms.
The target registry owns estimand-level theory. A study owns the treatment support
and whether outcome observation or an intermediate variable enters that theory. A
fresh record combines the two so identifying one study cannot mutate another target's
shared registry entry.
Every fact this needs is read off ``functional``, which
:func:`_point_functional` has already resolved. Reading the estimand and the data a
second time gave two implementations of the reference arm, and they printed different
answers into the same record.
Parameters
----------
base : Identification
The registered theory for the target.
design : PointTreatment
The declared column roles and mechanisms.
functional : BackdoorMeanContrast
The design-bound functional the study identified.
Returns
-------
Identification
The registered theory, narrowed to this design.
"""
target = functional.target
axis = functional.axis
missingness = functional.missingness
declared_level = functional.intermediate
direct = declared_level is not None
level = float(declared_level) if declared_level is not None else None
resolved: tuple[_ArmScope, str, str] | None = None
def arms() -> tuple[_ArmScope, str, str]:
"""Resolve the arm scope at most once, and only for a target that needs it."""
nonlocal resolved
if resolved is None:
resolved = _arm_scope(design, functional)
return resolved
assumptions: list[str] = []
for item in base.assumptions:
if item.startswith(MISSINGNESS_CAVEAT_PREFIX):
# The registry writes this as a condition because it does not know whether a
# response mechanism is declared. Drop it only where the design-bound MAR and
# response-positivity statements below replace it; otherwise the reader keeps
# the registry's own paragraph rather than nothing.
if design.missingness is None:
assumptions.append(item)
continue
if item.startswith(INTERMEDIATE_CAVEAT_PREFIX):
if design.intermediate is None:
assumptions.append(item)
continue
if target == "ey_obs":
if missingness is None and item.startswith(COMPLETE_OUTCOME_PREFIX):
assumptions.append(item)
continue
if target in {"att", "atc"} and item.startswith(REFERENCE_ONLY_POSITIVITY_PREFIX):
continue
if item.startswith(CONDITIONING_EVENT_PREFIX):
scope, _, _ = arms()
assumptions.append(_conditioning_event_statement(design, scope, target))
continue
if direct and item.startswith(CONSISTENCY_PREFIX):
assert level is not None
assumptions.append(
f"consistency: {design.outcome} = {design.outcome}(a, {level:g}) when "
f"{design.treatment_name} = a and {design.intermediate} = {level:g}"
)
continue
if direct and item.startswith(NO_CONFOUNDING_PREFIX):
assumptions.extend(
(
"no unmeasured treatment confounding: Y(a, z) is independent of A given W",
"no unmeasured intermediate confounding: Y(a, z) is independent of Z "
"given (A, W)",
)
)
continue
if item.startswith(TWO_ARM_POSITIVITY_PREFIX):
_, positivity, cells = arms()
assumptions.append(positivity)
if direct:
assert level is not None
assumptions.append(
f"intermediate positivity at {design.intermediate} = {level:g}: "
f"P({design.intermediate} = {level:g} | {design.treatment_name} = b, W) > 0 "
f"{cells}"
)
continue
assumptions.append(item)
nuisances = (
[OUTCOME_REGRESSION]
if target == "ey_obs" and missingness is not None
else []
if target == "ey_obs"
else list(base.required_nuisances)
)
if direct:
nuisances.append(INTERMEDIATE_MECHANISM)
if missingness is not None:
conditioning = (
f"({design.treatment_name}, {design.intermediate}, W)"
if direct
else f"({design.treatment_name}, W)"
)
assumptions.append(
f"missingness at random for {missingness}: {design.outcome} is independent "
f"of {missingness} given {conditioning}"
)
if target == "ey_obs":
assumptions.append(
f"response positivity for {missingness}: P({missingness} = 1 | "
f"{design.treatment_name}, W) > 0 almost surely"
)
elif target in POPULATION_INTERVENTION_TARGETS:
reference = functional.reference_arm
assumptions.append(
f"response positivity for {missingness}: P({missingness} = 1 | "
f"{design.treatment_name}, W) > 0 almost surely, for the natural-course term"
)
assumptions.append(
f"product positivity at the reference arm: P({design.treatment_name} = "
f"{reference!r} | W) P({missingness} = 1 | {design.treatment_name} = "
f"{reference!r}, W) > 0 almost surely"
)
elif axis == "policy":
assumptions.append(
f"response positivity for {missingness}: P({missingness} = 1 | "
f"{design.treatment_name} = d(a, W), W) > 0 almost surely wherever a declared "
"shift sends an observed dose"
)
elif direct or target in {"att", "atc", "ey", "ey1", "ey0"}:
_, _, cells = arms()
assumptions.append(
f"response positivity for {missingness}: P({missingness} = 1 | "
f"{design.treatment_name} = b, W) > 0 {cells}"
)
else:
assumptions.append(
f"response positivity for {missingness}: P({missingness} = 1 | "
f"{design.treatment_name} = a, W) > 0 almost surely wherever the target "
"functional evaluates arm a"
)
if direct:
assumptions.append(
f"observation-mechanism restriction: {missingness} is independent of "
f"{design.intermediate} given ({design.treatment_name}, W), so the fitted "
f"response mechanism excludes {design.intermediate}"
)
nuisances.append(MISSINGNESS_MECHANISM)
treatment_missingness = functional.treatment_missingness
if treatment_missingness is not None:
# The composite indicator's conditions: Y(a) independent of the treatment's
# recording given (A, W), and the composite probability bounded away from zero.
recorded = f"{treatment_missingness} = 1"
response = f", {missingness} = 1" if missingness is not None else ""
assumptions.append(
f"treatment missing at random for {treatment_missingness}: "
f"{design.outcome}(a) is independent of {treatment_missingness} given "
f"({design.treatment_name}, W); the recording may depend on "
f"{design.treatment_name}, so "
f"P({design.treatment_name} = a | W) is not identified and only "
f"P({design.treatment_name} = a | {recorded}, W) is"
)
assumptions.append(
f"composite positivity: P({design.treatment_name} = a, {recorded}{response} | W) > 0 "
"almost surely for every arm a"
)
nuisances.append(TREATMENT_OBSERVATION_MECHANISM)
# Definedness, stated for each estimand a boundary value leaves undefined. The three
# are one class of condition, and `cleverly.inference.delta` raises on each of them.
if target == "paf":
assumptions.append(
f"the natural-course risk E({design.outcome}) is strictly positive, so the "
"attributable fraction is defined"
)
if target == "rr":
assumptions.append(
f"both counterfactual risks are strictly positive: "
f"E_W[E({design.outcome} | {design.treatment_name} = a, W)] > 0 for every "
f"comparison arm a and at {design.treatment_name} = "
f"{functional.reference_arm!r}, so the log risk ratio is defined"
)
if target == "or":
assumptions.append(
f"both counterfactual risks lie strictly inside (0, 1): "
f"0 < E_W[E({design.outcome} | {design.treatment_name} = a, W)] < 1 for every "
f"comparison arm a and at {design.treatment_name} = "
f"{functional.reference_arm!r}, so the odds are finite"
)
if target == "ey_obs" and missingness is None:
assumptions.append(
"the natural-course mean is a functional of the observed law alone: neither "
f"the adjustment set {list(design.adjustment)} nor a treatment mechanism "
"enters it, which is why this record requires no nuisance estimate"
)
elif target == "ey_obs":
assumptions.append(
"the natural-course mean averages m(A, W) over the empirical joint law of "
"(A, W); the response and outcome regressions enter, but no treatment "
"mechanism or treatment-exchangeability assumption does"
)
references = base.references
if target == "ey_obs" and missingness is not None:
references = (*references, "Díaz, Carone & van der Laan (2016)")
if target in POPULATION_INTERVENTION_TARGETS and missingness is not None:
references = (
*references,
"Díaz, Carone & van der Laan (2016)",
"Díaz & van der Laan (2017)",
)
restrictions: list[str] = []
if design.treatment_adjustment is not None and set(design.treatment_adjustment) != set(
design.adjustment
):
restrictions.append(treatment_mechanism_restriction(design.treatment_adjustment))
return Identification(
assumptions=tuple(assumptions),
required_nuisances=tuple(nuisances),
dr_condition=_bound_dr_condition(base, design, functional, direct=direct),
references=references,
interpretations=base.interpretations,
model_restrictions=tuple(restrictions),
)
def _declared_interventions(actual: Any) -> tuple[str, InterventionKind] | None:
"""The field that holds a typed estimand's interventions, and their kind."""
if isinstance(actual, (RegimeMean, RegimeContrast)):
return "regimens", "regime"
if isinstance(actual, (ModifiedTreatmentPolicy, ModifiedTreatmentPolicyEffect)):
return "policies", "policy"
if isinstance(actual, (IncrementalMean, IncrementalEffect)):
return "interventions", "incremental"
return None
def _method_refusal(method: str, estimand: str, declared: Refusal | None) -> Refusal:
"""The refusal of ``estimate(method=...)``: the method record's refusal, named in context.
Parameters
----------
method : str
The requested method name.
estimand : str
The class name of the identified estimand.
declared : Refusal or None
The refusal of the method's :class:`~cleverly.methods.MethodAvailability` record.
``None`` means the study declares no record for the method.
Returns
-------
Refusal
The record's code, alternative and roadmap id, with a condition that names the method
and the estimand.
"""
if declared is None:
return Refusal(
RefusalReason.IMPLEMENTATION_GAP,
f"method {method!r} cannot estimate {estimand}: the method does not declare support.",
alternative="Choose a method that available_methods() lists as available.",
roadmap="X64",
)
return replace(
declared, condition=f"method {method!r} cannot estimate {estimand}: {declared.condition}"
)
_CONTINUOUS_DOSE_ALTERNATIVE = (
"Each of these pathwise-differentiable targets answers a different question. "
"ModifiedTreatmentPolicy([Shift(0.5, cap=None)]) is the mean outcome when every received "
"dose rises by 0.5. MSMProjection(MSM.linear(doses=[0.0, 0.5, 1.0])) projects the "
"dose-response curve at those doses onto a line."
)
def _continuous_dose_refusal(actual: Any, axis: str) -> Refusal:
"""The refusal of an estimand a continuous dose cannot index, by the estimand's axis.
Parameters
----------
actual : Any
The typed estimand.
axis : str
The estimand's parameter axis.
Returns
-------
Refusal
One record: the mean under a fixed dose or a deterministic rule is not pathwise
differentiable, a known stochastic density has no estimator here, and an odds tilt,
a risk-ratio tilt or another axis cannot index a dose at all.
"""
name = type(actual).__name__
declared = _declared_interventions(actual)
items = () if declared is None else tuple(getattr(actual, declared[0]))
if axis == "arm":
return Refusal(
RefusalReason.DIFFERENT_ESTIMAND,
f"{name} is indexed by treatment arm, and on a continuous dose an arm is one "
"fixed dose. The mean outcome under a fixed dose is not pathwise differentiable, "
"so no regular estimator of it exists.",
alternative=_CONTINUOUS_DOSE_ALTERNATIVE,
)
if axis == "regime" and any(isinstance(item, Stochastic) for item in items):
return Refusal(
RefusalReason.IMPLEMENTATION_GAP,
f"{name} sets the continuous dose by a known stochastic density. That target is "
"well posed, and this package has no estimator for a stochastic intervention on a "
"continuous dose.",
alternative=(
"A randomized modified policy with a finite known randomizer is supported: "
"ModifiedTreatmentPolicy([ModifiedPolicy(name, pieces={...}, "
"randomizer=Randomizer(values, probabilities), policy_kind='known')]). "
"MSMProjection(MSM.linear(doses=[0.0, 0.5, 1.0])) is the other route."
),
roadmap="X64",
)
if axis == "regime":
return Refusal(
RefusalReason.DIFFERENT_ESTIMAND,
f"{name} sets the continuous dose by a fixed value or a deterministic rule. That "
"parameter is not pathwise differentiable, so no regular estimator of it exists.",
alternative=_CONTINUOUS_DOSE_ALTERNATIVE,
)
if axis == "ipsi":
return Refusal(
RefusalReason.INVALID_DECLARATION,
f"{name} tilts the odds of treatment, and a continuous dose has no odds. Declare "
"a modified treatment policy.",
)
if axis == "rr_tilt":
return Refusal(
RefusalReason.INVALID_DECLARATION,
f"{name} is a risk-ratio tilt, which sets a binary treatment, and a continuous "
"dose is not binary.",
)
if axis == "learned_rule":
return Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"LearnedRuleValue over a continuous treatment is not implemented.",
alternative=(
"A fixed policy is available: ModifiedTreatmentPolicy([ModifiedPolicy(..., "
"policy_kind='known')])."
),
roadmap="F27",
)
return Refusal(
RefusalReason.INVALID_DECLARATION,
f"{name} is indexed by {axis}, which a continuous dose does not provide; a "
"continuous-dose design supports modified treatment policies and MSM projections.",
)
def _variant_roadmap(method: str, functional: Any, design: Any) -> str:
"""The roadmap item that holds a DR-TMLE or C-TMLE fit of a functional.
Parameters
----------
method : {"drtmle", "collaborative_tmle"}
The refused variant.
functional : Any
The identified functional, with its ``axis``, ``target`` and ``longitudinal``.
design : Any
The declared design, or ``None`` for an effect built without a study.
Returns
-------
str
F5 for the compositions it keeps, X50 or X51 for the ones that moved there, and X62
for the natural course with complete outcomes.
"""
axis, target = functional.axis, functional.target
# F5 keeps a longitudinal functional, a continuous dose, and ATT, ATC, incremental,
# modified-policy and risk-ratio-tilt fits under both variants.
if (
functional.longitudinal
or getattr(design, "treatment_kind", None) == "continuous"
or axis in {"ipsi", "policy", "rr_tilt"}
or (axis == "arm" and target in {"att", "atc"})
):
return "F5"
moved = {"drtmle": "X50", "collaborative_tmle": "X51"}[method]
if axis in {"regime", "msm"}:
return moved
natural_course = target == NATURAL_COURSE_TARGET or target in POPULATION_INTERVENTION_TARGETS
if natural_course and getattr(design, "missingness", None) is not None:
return moved
# X62 part (d): the natural course with complete outcomes.
return "X62"
def _point_functional(
design: PointTreatment,
data: CausalData,
estimand: PointEstimand,
) -> BackdoorMeanContrast:
"""Construct the canonical design-bound functional for one point estimand."""
actual, target = _point_target(estimand, data)
registered = TARGETS.get(target)
if registered is None:
raise CapabilityError(
f"{type(estimand).__name__} is not an evidenced point estimand",
reason=RefusalReason.INVALID_DECLARATION,
)
declared = _declared_interventions(actual)
interventions: tuple[Any, ...] = () if declared is None else tuple(getattr(actual, declared[0]))
if isinstance(actual, LearnedRuleValue):
# The rule rides where a regime would, so the engine adapter reads it off the
# functional as it reads every other axis's declaration.
interventions = (actual.rule,)
msm = actual.model if isinstance(actual, MSMProjection) else None
return BackdoorMeanContrast(
outcome=design.outcome,
treatment=design.treatment_name,
adjustment=tuple(design.adjustment),
target=target,
axis=registered.parameter_axis,
reference=getattr(actual, "reference", None),
interventions=interventions,
msm=msm,
intermediate=(
estimand.intermediate if isinstance(estimand, ControlledDirectEffect) else None
),
missingness=design.missingness,
treatment_missingness=design.treatment_missingness,
intermediate_name=design.intermediate,
treatment_levels=(tuple(data.treatment_levels) if _design_has_arms(design) else ()),
treatment_value=(actual.treatment if isinstance(actual, CounterfactualMean) else None),
)
def _matches_registered_point_identification(
actual: Identification,
registered: Identification,
identified: Any,
data: CausalData,
axis: str,
) -> bool:
"""Match a complete record against the one the registered design reconstructs."""
functional = getattr(identified, "functional", None)
study = getattr(identified, "_study", None)
design = getattr(study, "design", None)
if type(functional) is not BackdoorMeanContrast or type(design) is not PointTreatment:
return False
# Both calls sit inside the guard. A tampered estimand reaches the narrowing as well
# as the functional -- a controlled direct effect whose intermediate is not a number
# raises from ``float(...)`` -- and an exception out of a provenance matcher is a
# different failure from the refusal it is asked for.
try:
expected_functional = _point_functional(design, data, identified.estimand)
expected_registered = TARGETS[expected_functional.target]
if registered != expected_registered.identification or axis != expected_functional.axis:
return False
expected = _point_identification(registered, design, expected_functional)
except (AttributeError, CapabilityError, DataError, KeyError, TypeError, ValueError):
return False
# ``actual == expected`` compares every field of the record, ``model_restrictions``
# included, so a design-level restriction dropped or added after the fact is refused.
return functional == expected_functional and actual == expected
_LONGITUDINAL_IDENTIFICATION = Identification(
assumptions=(
"consistency: each observed history equals the potential history under its "
"realized regimen",
"no interference: one unit's potential history does not depend on other units' regimens",
"sequential exchangeability given the recorded history at every node",
"sequential positivity for treatment and remaining under observation",
),
required_nuisances=("sequential_outcome_regressions", "treatment_and_censoring_mechanisms"),
dr_condition=(
"the longitudinal remainder is a sum of stagewise products of outcome-regression "
"and cumulative mechanism errors"
),
references=("van der Laan & Gruber (2012)",),
)
#: The identification of a regimen set with a known policy node. Positivity is required
#: wherever the plan puts mass, which at a policy node is every arm the policy can draw,
#: and the policy history excludes the current arm, so Assumption 2 of Díaz, Williams,
#: Hoffman and Schenck (2023) suffices.
_LONGITUDINAL_POLICY_IDENTIFICATION = Identification(
assumptions=(
*_LONGITUDINAL_IDENTIFICATION.assumptions[:3],
"sequential positivity for treatment and remaining under observation, wherever a "
"regimen or its known policy puts mass",
"each policy density is known and fixed before the fit, and reads the history "
"before its node but not the treatment at that node",
),
required_nuisances=_LONGITUDINAL_IDENTIFICATION.required_nuisances,
dr_condition=_LONGITUDINAL_IDENTIFICATION.dr_condition,
references=(
"van der Laan & Gruber (2012)",
"Díaz, Williams, Hoffman & Schenck (2023)",
),
)
#: The two readings of the longitudinal policy functional (X40 part (d)). Díaz, Williams,
#: Hoffman and Schenck (2023), Theorem 1: under strong sequential randomization (Assumption 3)
#: the functional identifies the LMTP, and under standard sequential randomization
#: (Assumption 2) it identifies the LMTP-SI. Williams and Díaz (2025) correct the software
#: paper's exchangeability to the strong form.
_LONGITUDINAL_MTP_INTERPRETATIONS = (
Interpretation(
label="modified treatment policy effect",
meaning=(
"the mean outcome had each node's natural dose a_t been replaced by "
"d_t(a_t, h_t), with the history the earlier replacements produce"
),
assumptions=(
"strong sequential randomization: at every node, every common cause of the dose, "
"the outcome, and the later covariates and doses is in the recorded history "
"(Díaz et al. 2023, Assumption 3; the strong exchangeability of Williams & Díaz "
"2025, A_t independent of Y(a_t, d_{t+1}(A_{t+1}(a_t), H_{t+1}(a_t))) given H_t)",
),
references=(
"Díaz, Williams, Hoffman & Schenck (2023), Theorem 1 (ii) and Assumption 3",
"Williams & Díaz (2025)",
),
),
Interpretation(
label="stochastic intervention (LMTP-SI)",
meaning=(
"the mean outcome had each node's dose been drawn from the law of d_t(A_t, H_t) "
"given the history"
),
assumptions=(
"standard sequential randomization: at every node, every common cause of the dose "
"and the next covariates, which at the last node are the outcome, is in the "
"recorded history (Díaz et al. 2023, Assumption 2)",
),
references=("Díaz, Williams, Hoffman & Schenck (2023), Theorem 1 (i) and Assumption 2",),
),
)
#: The identification of a regimen set with a modified treatment policy node. The shared
#: assumptions are what both readings of :data:`_LONGITUDINAL_MTP_INTERPRETATIONS` need,
#: and one line that names the exchangeability each reading adds.
_LONGITUDINAL_MTP_IDENTIFICATION = Identification(
assumptions=(
*_LONGITUDINAL_IDENTIFICATION.assumptions[:2],
"sequential exchangeability in the form the chosen interpretation names: strong "
"sequential randomization (Díaz et al. 2023, Assumption 3) for the modified "
"treatment policy effect, or standard sequential randomization (Assumption 2) for "
"the LMTP-SI of their Theorem 1 (i)",
"positivity: the dose a policy assigns, d(a_t, h_t), lies in the support of A_t given "
"h_t wherever a_t does (Díaz et al. 2023, Assumption 1, which Williams & Díaz 2025 "
"leave unchanged), and so does remaining under observation",
"each policy is known and fixed before the fit, or its cap is a declared tuner with a "
"declared target (X40 part (e)), and on a continuous node it is piecewise smooth "
"invertible (Assumption 4); a policy reads the history and the unit's own treatment "
"at its node, and its earlier treatment columns stand for the intervened values",
),
required_nuisances=_LONGITUDINAL_IDENTIFICATION.required_nuisances,
dr_condition=(
"the remainder is a sum over nodes of the products of the density ratio error and "
"the outcome regression error (Díaz et al. 2023, Theorem 3), with bounded ratios"
),
references=("Díaz, Williams, Hoffman & Schenck (2023)", "Williams & Díaz (2025)"),
interpretations=_LONGITUDINAL_MTP_INTERPRETATIONS,
)
#: The identification of a regimen set with an odds tilt node. Kennedy (2019), Theorem 1:
#: under consistency and exchangeability (Assumptions 1 and 2) the tilt is identified with
#: no treatment positivity at its nodes. Under dropout, Kim, Kennedy and Naimi (2021) add
#: missingness at random (Assumption A2-M) and positivity of remaining under observation
#: (Assumption A3). A label, rule or policy node beside a tilt keeps its own positivity.
_LONGITUDINAL_TILT_IDENTIFICATION = Identification(
assumptions=(
*_LONGITUDINAL_IDENTIFICATION.assumptions[:3],
"no treatment positivity at a tilt node: the tilt keeps each arm's support "
"(Kennedy 2019, Theorem 1). A label, rule or policy node needs positivity wherever "
"its plan puts mass",
"dropout is missing at random given the recorded history, and each unit at risk "
"remains under observation with positive probability (Kim, Kennedy & Naimi 2021, "
"Assumptions A2-M and A3)",
),
required_nuisances=_LONGITUDINAL_IDENTIFICATION.required_nuisances,
dr_condition=(
"not doubly robust: at each tilt node the remainder holds the product of the "
"treatment-mechanism error with the sum of the outcome-regression and "
"treatment-mechanism errors (Kennedy 2019, Section 3.3 and Theorem 3)"
),
references=("Kennedy (2019)", "Kim, Kennedy & Naimi (2021)"),
)
[docs]
@dataclass(frozen=True)
class ExplicitAdjustmentProvider:
"""Identify effects from declared adjustment sets or sequential histories.
Parameters
----------
name : str
Provider name recorded in the identified effect.
"""
name: str = "explicit-adjustment"
[docs]
def identify(self, study: CausalStudy, estimand: PointEstimand) -> IdentifiedEffect:
"""Bind a supported estimand to the design's observed-data functional.
Parameters
----------
study : CausalStudy
Validated observed data and design.
estimand : Estimand
Typed causal quantity.
Returns
-------
IdentifiedEffect
Functional, assumptions, and method availability for the question.
"""
design = study.design
if isinstance(design, (LongitudinalTreatment, TimeToEvent, LongPanel)):
return self._identify_longitudinal(study, estimand)
return self._identify_point(study, estimand)
def _identify_point(self, study: CausalStudy, estimand: PointEstimand) -> IdentifiedEffect:
design = study.design
assert isinstance(design, PointTreatment)
data = study.data
assert isinstance(data, CausalData)
actual, target = _point_target(estimand, data)
if (
isinstance(actual, CounterfactualMean)
and actual.treatment is not None
and actual.treatment not in data.treatment_levels
):
raise DataError(
f"treatment {actual.treatment!r} is not available; "
f"choose from {list(data.treatment_levels)}"
)
if target not in TARGETS:
raise CapabilityError(
f"{type(estimand).__name__} is not an evidenced point estimand",
reason=RefusalReason.INVALID_DECLARATION,
)
if isinstance(actual, LearnedRuleValue):
# The rows of the learned-rule refusal table that the data decide, before the generic
# axis and design checks below: their sentences name other remedies.
refuse_learned_rule_composition(data)
# The kinds of the items are a property of the estimand, so they are checked before
# the design. Not in ``_point_functional``: the provenance matcher calls that inside
# a ``try`` that reads a refusal as a mismatch.
declared = _declared_interventions(actual)
if declared is not None:
field_name, kind = declared
items = getattr(actual, field_name)
holder = f"{type(actual).__name__}.{field_name}"
# Refused rather than read: an iterator is empty after its first read, and
# ``_point_functional`` and the provenance matcher each read the set again.
problem = None
if isinstance(items, Mapping):
problem = (
"It is a mapping, and a point-treatment design reads no mapping key. Give "
"each item its name with name="
)
elif isinstance(items, Iterator):
problem = (
f"It is an iterator, {type(items).__name__}, and an iterator is empty "
"after its first read"
)
elif (
isinstance(items, (str, bytes))
or getattr(items, "ndim", 1) == 0
or not isinstance(items, Iterable)
):
problem = f"It received one item, {items!r}. Write one item as a one-item tuple"
if problem is not None:
raise DataError(f"{holder} must be a sequence, such as a tuple. {problem}")
refuse_mixed_interventions(items, kind=kind, holder=holder)
if data.has_missing_treatment:
refusal = missing_treatment_design_refusal(
data,
target=target,
axis=TARGETS[target].parameter_axis,
intermediate=design.intermediate is not None,
)
if refusal is not None:
raise refusal
axis = TARGETS[target].parameter_axis
# The rule is about the parameter *axis*, not about a list of target names. Naming
# targets refused `msm` -- whose axis is `msm`, indexed by working-model term, and
# which the density-ratio path fits happily on a dose (tests/unit/test_continuous_msm.py)
# -- while telling the caller it was "arm-indexed", which it is not. A regime is a
# density over arms and an incremental intervention tilts an odds, so those two stay
# refused; both statements now come from the axis rather than from a spelling.
if design.support == "vector" and axis != "policy":
reason = (
"MSMProjection projects the dose-response curve of one dose over a grid on "
"that dose, and a vector has no one dose axis and no integration grid here"
if axis == "msm"
else f"it is indexed by {axis}, and such a treatment has no arms"
)
raise CapabilityError(
f"treatment {design.treatment_name!r} is a vector with a continuous "
f"component, and {type(actual).__name__} is refused: {reason}.",
reason=(
RefusalReason.DIFFERENT_ESTIMAND
if axis in {"arm", "regime"}
else RefusalReason.INVALID_DECLARATION
),
alternative=(
"A vector design supports ModifiedTreatmentPolicy and "
"ModifiedTreatmentPolicyEffect with policies=[JointPolicy(...), ...]."
),
)
if design.treatment_kind == "continuous" and axis not in {"policy", "msm"}:
raise CapabilityError.from_refusal(_continuous_dose_refusal(actual, axis))
if getattr(actual, "horizons", None) is not None:
raise CapabilityError(
f"{type(actual).__name__}(horizons=...) selects follow-up times from a "
"sequential fit, and a point-treatment design has one time point; declare "
"LongitudinalTreatment to ask that question",
reason=RefusalReason.INVALID_DECLARATION,
)
if isinstance(actual, MSMProjection) and actual.regimens is not None:
raise CapabilityError(
"MSMProjection(regimens=...) projects over longitudinal regimen cells; on a "
"point-treatment design the projection is over the declared arms or doses. "
"Drop regimens=, or declare LongitudinalTreatment",
reason=RefusalReason.INVALID_DECLARATION,
)
reference = getattr(actual, "reference", None)
if (
reference is not None
and _design_has_arms(design)
and axis == "arm"
and reference not in data.treatment_levels
):
raise DataError(
f"reference {reference!r} is not a treatment level; "
f"available: {list(data.treatment_levels)}"
)
if isinstance(estimand, ControlledDirectEffect):
if design.intermediate is None:
raise CapabilityError(
"ControlledDirectEffect needs PointTreatment(intermediate=...)",
reason=RefusalReason.INVALID_DECLARATION,
)
if estimand.intermediate not in (0.0, 1.0):
raise DataError("the evidenced binary intermediate must be fixed at 0 or 1")
elif design.intermediate is not None:
raise CapabilityError(
"a design with intermediate= must identify ControlledDirectEffect explicitly",
reason=RefusalReason.INVALID_DECLARATION,
)
if isinstance(actual, MSMProjection):
refuse_continuous_msm_mechanisms(
data, subject="MSMProjection", missingness="PointTreatment(missingness=...)"
)
# The estimand and the design decide these data conditions, so they refuse here.
check_stratified_targets(
data,
incremental=isinstance(actual, (IncrementalMean, IncrementalEffect)),
msm=actual.model if isinstance(actual, MSMProjection) else None,
)
functional = _point_functional(design, data, estimand)
return IdentifiedEffect(
estimand=estimand,
functional=functional,
identification=_point_identification(
TARGETS[target].identification, design, functional
),
provider=self,
_study=study,
)
def _identify_longitudinal(
self, study: CausalStudy, estimand: PointEstimand
) -> IdentifiedEffect:
design = study.design
assert isinstance(design, (LongitudinalTreatment, TimeToEvent, LongPanel))
if isinstance(estimand, LearnedRuleValue):
raise CapabilityError(
"LearnedRuleValue is a point-treatment target, and a learned longitudinal "
"regimen has no supported path. Van der Laan and Luedtke (2015) give one at "
"two time points, and X11 (b) in docs/roadmap.md holds it.",
reason=RefusalReason.IMPLEMENTATION_GAP,
alternative="DynamicRegimen(..., rule_kind='known') fits a fixed regimen.",
roadmap="X11",
)
if isinstance(estimand, (IncrementalMean, IncrementalEffect)):
regimens: Any = _longitudinal_tilts(estimand)
horizons = estimand.horizons
elif isinstance(estimand, (RegimeMean, RegimeContrast, MSMProjection)):
regimens = estimand.regimens
horizons = estimand.horizons
else:
raise CapabilityError(
f"{type(estimand).__name__} is not a longitudinal regimen functional",
reason=RefusalReason.INVALID_DECLARATION,
)
if regimens is None:
raise DataError("a longitudinal estimand must declare regimens")
target = "msm_regimen" if isinstance(estimand, MSMProjection) else estimand.name
data = study.data
assert isinstance(data, LongitudinalData)
# The conditioning set of the treatment mechanism at each decision, in the column
# order of ``LongitudinalData.history_design``: the covariate history up to that
# node, then every earlier treatment. A held design has one decision, at node 1.
history = tuple(
(node, (*data.history_names(time), *data.treatment_names[: time - 1]))
for time, node in enumerate(data.treatment_names, start=1)
)
functional = BackdoorMeanContrast(
outcome=design.outcome_role if isinstance(design, LongPanel) else design.outcome,
treatment=tuple(data.treatment_names),
adjustment=tuple(design.baseline),
target=target,
axis="regimen" if target != "msm_regimen" else "msm",
reference=getattr(estimand, "reference", None),
interventions=regimens,
horizons=None if horizons is None else tuple(horizons),
msm=estimand.model if isinstance(estimand, MSMProjection) else None,
longitudinal=True,
history=history,
)
return IdentifiedEffect(
estimand=estimand,
functional=functional,
identification=(
_LONGITUDINAL_TILT_IDENTIFICATION
if declares_tilt(regimens)
else _LONGITUDINAL_MTP_IDENTIFICATION
if declares_mtp(regimens)
else _LONGITUDINAL_POLICY_IDENTIFICATION
if declares_policy(regimens)
else _LONGITUDINAL_IDENTIFICATION
),
provider=self,
_study=study,
)
def _longitudinal_tilts(estimand: IncrementalMean | IncrementalEffect) -> dict[str, Incremental]:
"""The plans of an incremental estimand on a longitudinal design, one per tilt.
Each tilt sets every node, under its own ``name``, as ``IPSI([...])`` reads a sequence.
"""
holder = f"{type(estimand).__name__}.interventions"
items = estimand.interventions
if isinstance(items, (Mapping, Iterator, str, bytes)) or not isinstance(items, Iterable):
raise DataError(
f"{holder} must be a sequence of Incremental, such as a tuple. It received "
f"{type(items).__name__}"
)
items = tuple(items)
refuse_mixed_interventions(items, kind="incremental", holder=holder)
names = [item.name for item in items]
repeated = sorted({name for name in names if names.count(name) > 1})
if repeated:
raise DataError(
f"two tilts share the name {repeated[0]!r}; a name labels a parameter, so give "
"each Incremental its own name="
)
return {item.name: item for item in items}
[docs]
class CausalStudy:
"""Validate observed data against an immutable causal-study design.
Parameters
----------
data : dataframe, CausalData, or LongitudinalData
Observed study data. Pandas and Polars dataframes are supported.
design : PointTreatment, LongitudinalTreatment, TimeToEvent, or LongPanel
Column roles and treatment-time structure.
protocol : StudyProtocol or None
Scientific study protocol. ``None`` records that no protocol was supplied.
Attributes
----------
data : CausalData or LongitudinalData
design : PointTreatment, LongitudinalTreatment, TimeToEvent, or LongPanel
protocol : StudyProtocol or None
See Also
--------
PointTreatment : The design declaration for treatment given once.
LongitudinalTreatment : The design declaration for time-varying treatment.
TimeToEvent : The design declaration for one follow-up time and one event code per unit.
LongPanel : The design declaration for one row per unit and decision time.
IdentifiedEffect : What :meth:`identify` returns.
Examples
--------
>>> from cleverly import ATE, CausalStudy, PointTreatment
>>> from cleverly.datasets import make_linear_ate
>>> frame, _ = make_linear_ate(n=40, seed=1)
>>> study = CausalStudy(
... frame,
... design=PointTreatment(
... outcome="Y", treatment="A", adjustment=("W1", "W2", "W3", "W4")
... ),
... )
>>> study.identify(ATE()).estimand.name
'ate'
"""
_protocol: StudyProtocol | None
def __init__(
self,
data: Any,
*,
design: PointTreatment | LongitudinalTreatment | TimeToEvent | LongPanel,
protocol: StudyProtocol | None = None,
) -> None:
if protocol is not None and not isinstance(protocol, StudyProtocol):
raise TypeError("protocol must be a StudyProtocol or None")
self._design = design
self._data = design.prepare(data)
self._protocol = protocol
@property
def design(self) -> PointTreatment | LongitudinalTreatment | TimeToEvent | LongPanel:
"""Return the immutable study design."""
return self._design
@property
def data(self) -> CausalData | LongitudinalData:
"""Return the validated internal data representation."""
return self._data
@property
def protocol(self) -> StudyProtocol | None:
"""Return the scientific study protocol, if one was supplied."""
return self._protocol
[docs]
def identify(
self,
estimand: PointEstimand,
*,
provider: IdentificationProvider | None = None,
) -> IdentifiedEffect:
"""Bind one typed estimand to an observed-data functional.
Parameters
----------
estimand : Estimand
Typed causal quantity, such as ``ATE()``. String aliases are not accepted.
provider : IdentificationProvider or None
Identification strategy. ``None`` uses explicit adjustment or history.
Returns
-------
IdentifiedEffect
Question, functional, assumptions, and available methods before fitting.
Raises
------
CapabilityError
If the estimand type or its composition with the design or the data is
unsupported, or if a typed estimand holds an intervention of another kind.
DataError
If the estimand names a value that the data lack, or if the set of a point
estimand is a mapping, an iterator, a string, or a single item.
See Also
--------
IdentifiedEffect : The identified question returned by this method.
CausalStudy.estimate : Identify and estimate in one call.
Notes
-----
Pass a typed estimand such as ``ATE()``. String aliases are refused with a
:class:`CapabilityError` so the failure identifies the required public object.
Examples
--------
>>> from cleverly import ATE, CausalStudy, PointTreatment
>>> from cleverly.datasets import make_linear_ate
>>> frame, _ = make_linear_ate(n=40, seed=1)
>>> study = CausalStudy(
... frame,
... design=PointTreatment(
... outcome="Y", treatment="A", adjustment=("W1", "W2", "W3", "W4")
... ),
... )
>>> effect = study.identify(ATE(reference=0))
>>> effect.estimand.reference
0
"""
if not isinstance(estimand, PointEstimand):
raise CapabilityError(
f"{estimand!r} is not a typed causal estimand.",
reason=RefusalReason.API_REDIRECTION,
alternative=(
"Pass one of the objects exported by cleverly, such as ATE()"
+ (
f". For {estimand!r}, that is {_STRING_ESTIMANDS[estimand]}"
if isinstance(estimand, str) and estimand in _STRING_ESTIMANDS
else ""
)
+ "."
),
)
effect = (provider or ExplicitAdjustmentProvider()).identify(self, estimand)
return replace(effect, protocol=self.protocol, _study=self)
[docs]
def estimate(
self,
estimand: PointEstimand,
method: str | EstimationMethod = "tmle",
**overrides: Any,
) -> TMLEResult | LongitudinalResult:
"""Identify and estimate one causal quantity.
Parameters
----------
estimand : Estimand
Typed causal quantity.
method : str or EstimationMethod
Named method or configured method object.
**overrides : Any
Method fields overridden for this fit.
Returns
-------
TMLEResult or LongitudinalResult
Fitted result for the declared design.
"""
return self.identify(estimand).estimate(method=method, **overrides)
[docs]
@dataclass(frozen=True)
class IdentifiedEffect: # numpydoc ignore=PR01
"""Bind an estimand to a functional and its identification assumptions.
Parameters
----------
estimand : Estimand
Typed causal quantity.
functional : BackdoorMeanContrast
Normalized observed-data target.
identification : Identification
Assumptions, nuisance requirements, and remainder condition.
provider : IdentificationProvider
Provider that performed identification.
protocol : StudyProtocol or None
Scientific study protocol stamped by :class:`CausalStudy`.
See Also
--------
CausalStudy.identify : What returns this object.
MethodAvailability : One row of :meth:`available_methods`.
cleverly.EstimationMethod : The configuration :meth:`estimate` accepts.
Notes
-----
``_study`` is the bound study, kept private because it is not part of the reported
causal question. An effect built without a study leaves it ``None``, and such an
effect cannot be refitted.
Examples
--------
:meth:`CausalStudy.identify` returns this object. It carries the assumptions the
estimate will rest on, and it is what an estimation method is passed to.
>>> from cleverly import ATE, CausalStudy, PointTreatment
>>> from cleverly.datasets import make_linear_ate
>>> frame, _ = make_linear_ate(n=40, seed=1)
>>> study = CausalStudy(
... frame,
... design=PointTreatment(
... outcome="Y", treatment="A", adjustment=("W1", "W2", "W3", "W4")
... ),
... )
>>> effect = study.identify(ATE())
>>> effect.estimand.name
'ate'
>>> sorted(method.name for method in effect.available_methods())[:2]
['collaborative_tmle', 'drtmle']
"""
estimand: PointEstimand
functional: BackdoorMeanContrast
identification: Identification
provider: IdentificationProvider
protocol: StudyProtocol | None = None
_study: CausalStudy | None = field(repr=False, compare=False, default=None)
[docs]
def available_methods(self) -> tuple[MethodAvailability, ...]:
"""Report estimation methods available for this functional.
Returns
-------
tuple of MethodAvailability
One capability record for each named method, including refusal reasons.
Notes
-----
The result includes unavailable methods and their refusal reasons. Selecting one
raises :class:`CapabilityError` before nuisance fitting starts.
``available`` is per method name. A configuration of an available method can still
refuse before any learner. For :class:`LearnedRuleValue`, the default
``TMLEMethod()`` refuses (X11 (e) in ``docs/roadmap.md``), and
``TMLEMethod(cross_fitting=CrossFitting(fold_evaluation=True))`` fits.
"""
point = not self.functional.longitudinal
target = self.functional.target
blocker = (
Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"a controlled direct effect fixes an intermediate variable, and neither "
"variant's score is derived for that functional.",
alternative="Fit method='tmle'.",
roadmap="F5",
)
if self.functional.intermediate is not None
else Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"no collaborative or doubly robust learned-rule result was reviewed.",
alternative="Fit method='tmle'.",
roadmap="F27",
)
if self.functional.axis == "learned_rule"
else None
)
variants = (
blocker is None
and point
and self.functional.axis == "arm"
and target in {"ate", "ey", "ey1", "ey0", "rr", "or"}
)
# Read off the declared design rather than off a fitted result, so a caller who
# asks which methods are available before fitting is told the same thing the fit
# would tell them. An effect built without a study cannot answer this, and it
# claims nothing.
design = getattr(self._study, "design", None)
# A design that declares its mechanism leaves C-TMLE nothing to select. A design that
# declares clusters keeps the method: strategy="derivative" (C-TMLE0) selects
# nothing and has the ordinary TMLE's clustered result, and every other strategy
# refuses id= before its first learner.
known = getattr(design, "treatment_probabilities", None) is not None
collaborative_item = _variant_roadmap("collaborative_tmle", self.functional, design)
drtmle_item = _variant_roadmap("drtmle", self.functional, design)
collaborative_blocker = (
KNOWN_CTMLE_REFUSAL
if known
else blocker
or Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"no collaborative score is evidenced for this functional.",
alternative="Fit method='tmle'.",
roadmap=collaborative_item,
)
)
collaborative = variants and not known
drtmle_blocker = blocker or Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"no reduced-dimension correction is evidenced for this functional.",
alternative="Fit method='tmle'.",
roadmap=drtmle_item,
)
if isinstance(design, PointTreatment) and design.support == "vector":
# A vector with a dose has no arms, as a continuous dose has none; each variant
# refuses it with the reason its engine gives a continuous treatment.
from .estimators.ctmle import vector_discrete_treatment_refusal
from .estimators.drtmle import drtmle_axis_refusal
collaborative_blocker = vector_discrete_treatment_refusal(design.treatment_name)
drtmle_blocker = drtmle_axis_refusal("policies")
tmle, tmle_reason, sdr, sdr_reason, ipsi, ipsi_reason = self._sequential_methods()
gcomp, gcomp_reason, ipw, ipw_reason = self._baseline_methods()
gformula, gformula_reason = self._gformula_method()
return (
MethodAvailability("tmle", tmle, tmle_reason),
MethodAvailability(
"collaborative_tmle",
collaborative,
None if collaborative else collaborative_blocker,
),
MethodAvailability("drtmle", variants, None if variants else drtmle_blocker),
MethodAvailability("sdr", sdr, sdr_reason),
MethodAvailability("ipsi", ipsi, ipsi_reason),
MethodAvailability("gcomp", gcomp, gcomp_reason),
MethodAvailability("ipw", ipw, ipw_reason),
MethodAvailability("parametric_gformula", gformula, gformula_reason),
MethodAvailability(
"riesz_tmle",
False,
Refusal(
RefusalReason.IMPLEMENTATION_GAP,
"the direct-Riesz engine and intervention-state representer are not "
"implemented.",
alternative="Fit method='tmle'.",
roadmap="R1",
),
),
MethodAvailability(
"ep",
False,
Refusal(
RefusalReason.DIFFERENT_ESTIMAND,
"EP requires a conditional-contrast estimand and its distinct result family. "
"P1 in docs/roadmap.md proposes the EP learner and that estimand.",
alternative=(
"To estimate the ATE inside each level of a discrete baseline "
"modifier, declare PointTreatment(..., strata=['V']) and request ATE()."
),
),
),
)
def _sequential_methods(
self,
) -> tuple[bool, Refusal | None, bool, Refusal | None, bool, Refusal | None]:
"""Whether ``tmle``, ``sdr`` and ``ipsi`` fit this functional, each with its refusal.
An odds tilt reads the treatment mechanism. ``LTMLE`` fits it only as a known
policy, when every tilt node's mechanism is declared (R1), and ``SDR`` never fits it
(R2). ``IPSI`` fits every longitudinal plan set, with or without a tilt node. On a
repeated outcome, a method is available only when its kind is in the allowlist of
estimators with a per-horizon reduction (``REPEATED_OUTCOME_KINDS``).
"""
if not self.functional.longitudinal:
return True, None, False, SDR_POINT_REFUSAL, False, IPSI_POINT_REFUSAL
data = getattr(self._study, "data", None)
tmle, tmle_reason, sdr, sdr_reason, ipsi, ipsi_reason = self._tilt_methods(data)
# The fit refuses a repeated outcome on a kind outside the allowlist before any
# learner (R7), so the listing must not offer it, and its reason is the fit's.
if getattr(data, "is_repeated", False):
kinds = _sequential.REPEATED_OUTCOME_KINDS
if "tmle" not in kinds:
tmle, tmle_reason = False, _sequential.repeated_outcome_refusal("LTMLE")
if "sdr" not in kinds:
sdr, sdr_reason = False, _sequential.repeated_outcome_refusal("SDR")
if "ipsi" not in kinds:
ipsi, ipsi_reason = False, _sequential.repeated_outcome_refusal("IPSI")
return tmle, tmle_reason, sdr, sdr_reason, ipsi, ipsi_reason
def _tilt_methods(
self, data: Any
) -> tuple[bool, Refusal | None, bool, Refusal | None, bool, Refusal | None]:
"""The three sequential methods' availability as the plan set's odds tilts allow."""
regimens = self.functional.interventions
# IPSI estimates no density ratio, so a plan with a policy node has no IPSI fit.
ipsi, ipsi_reason = (False, IPSI_POLICY_REFUSAL) if declares_mtp(regimens) else (True, None)
if not declares_tilt(regimens):
return True, None, True, None, ipsi, ipsi_reason
known = getattr(data, "known_mechanisms", None)
n_times = getattr(data, "n_times", 0)
declared = known is not None and all(
known.treatment_at(time) is not None for time in tilt_nodes(regimens, n_times)
)
return (
declared,
None if declared else LTMLE_TILT_REFUSAL,
False,
SDR_TILT_REFUSAL,
ipsi,
ipsi_reason,
)
def _baseline_methods(self) -> tuple[bool, Refusal | None, bool, Refusal | None]:
"""Whether ``gcomp`` and ``ipw`` fit this functional, each with its reason.
Both baselines fit longitudinal regimen means and contrasts only. They refuse a
working model over regimens and an odds tilt whose mechanism is estimated. ``ipw``
also refuses a modified treatment policy, which needs a density ratio it does not fit.
On a repeated outcome each is available only when its kind is in
``REPEATED_OUTCOME_KINDS``, as the fit refuses it before any learner.
"""
if not self.functional.longitudinal:
return False, BASELINE_POINT_REFUSAL, False, BASELINE_POINT_REFUSAL
if self.functional.msm is not None:
return False, BASELINE_MSM_REFUSAL, False, BASELINE_MSM_REFUSAL
gcomp, gcomp_reason, ipw, ipw_reason = self._baseline_plan_methods()
if getattr(getattr(self._study, "data", None), "is_repeated", False):
kinds = _sequential.REPEATED_OUTCOME_KINDS
if "gcomp" not in kinds:
gcomp = False
gcomp_reason = _sequential.repeated_outcome_refusal("SequentialGComputation")
if "ipw" not in kinds:
ipw, ipw_reason = False, _sequential.repeated_outcome_refusal("SequentialIPW")
return gcomp, gcomp_reason, ipw, ipw_reason
def _gformula_method(self) -> tuple[bool, Refusal | None]:
"""Whether ``parametric_gformula`` fits this functional, with its reason.
It fits longitudinal regimen means and contrasts only, and refuses a working model over
regimens. It accepts every plan kind, an odds tilt included, because it models the
natural treatment itself. On a repeated outcome it is unavailable, as the fit refuses
it before any learner.
"""
if not self.functional.longitudinal:
return False, BASELINE_POINT_REFUSAL
if self.functional.msm is not None:
return False, BASELINE_MSM_REFUSAL
if getattr(getattr(self._study, "data", None), "is_repeated", False):
kinds = _sequential.REPEATED_OUTCOME_KINDS
if "parametric_gformula" not in kinds:
return False, _sequential.repeated_outcome_refusal("ParametricGFormula")
return True, None
def _baseline_plan_methods(self) -> tuple[bool, Refusal | None, bool, Refusal | None]:
"""The two baselines' availability as the plan set allows."""
regimens = self.functional.interventions
if declares_tilt(regimens):
data = getattr(self._study, "data", None)
known = getattr(data, "known_mechanisms", None)
n_times = getattr(data, "n_times", 0)
declared = known is not None and all(
known.treatment_at(time) is not None for time in tilt_nodes(regimens, n_times)
)
if not declared:
return False, BASELINE_TILT_REFUSAL, False, BASELINE_TILT_REFUSAL
if declares_mtp(regimens):
return True, None, False, IPW_MTP_REFUSAL
return True, None, True, None
[docs]
def summary(self, *, protocol: ProtocolDetail = "full") -> str:
"""Return a readable identification summary.
Parameters
----------
protocol : {"full", "fingerprint"}
``"full"`` prints every field of the protocol record. ``"fingerprint"`` prints
only its schema and fingerprint, for a reader who has the record already.
Returns
-------
str
A printable block naming the estimand, the functional, and the assumptions.
"""
assumptions = "\n".join(f" - {item}" for item in self.identification.assumptions)
readings = [
line
for reading in self.identification.interpretations
for line in (
f" - {reading.label}: {reading.meaning}",
*(f" needs {item}" for item in reading.assumptions),
*(f" source: {item}" for item in reading.references),
)
]
restrictions = (
[
"model restrictions:",
*(f" - {item}" for item in self.identification.model_restrictions),
]
if self.identification.model_restrictions
else []
)
history = self.functional.history
# A longitudinal functional conditions each node's mechanism on its own history,
# so one list of the baseline covariates would understate every later node.
adjustment = (
[f"history at {node}: {list(names)}" for node, names in history]
if history
else [f"adjustment/history: {list(self.functional.adjustment)}"]
)
return "\n".join(
[
self.estimand.definition,
f"identified by {self.provider.name}: {self.functional.expression}",
*adjustment,
f"required nuisances: {list(self.identification.required_nuisances)}",
f"assumptions:\n{assumptions}",
*(
["interpretations, each with the assumption it adds:", *readings]
if readings
else []
),
*restrictions,
*protocol_summary_lines(self.protocol, protocol),
]
)
[docs]
def summary_lines(self, *, protocol: ProtocolDetail = "full") -> tuple[str, ...]:
"""Return identification facts for inclusion in result summaries.
Parameters
----------
protocol : {"full", "fingerprint"}
``"full"`` gives every field of the protocol record. ``"fingerprint"`` gives
only its schema and fingerprint.
Returns
-------
tuple of str
The same facts as lines, for a result summary to append.
"""
return (
f"causal estimand: {self.estimand.definition}",
f"identification: {self.provider.name}; {self.functional.expression}",
"required nuisances: " + ", ".join(self.identification.required_nuisances),
"identification assumptions: " + "; ".join(self.identification.assumptions),
*(
f"interpretation {reading.label}: {reading.meaning}; needs "
+ "; ".join(reading.assumptions)
for reading in self.identification.interpretations
),
*(
("model restrictions: " + "; ".join(self.identification.model_restrictions),)
if self.identification.model_restrictions
else ()
),
*protocol_summary_lines(self.protocol, protocol),
)
[docs]
def estimate(
self,
method: str | EstimationMethod = "tmle",
**overrides: Any,
) -> TMLEResult | LongitudinalResult:
"""Estimate the identified functional.
Parameters
----------
method : str or EstimationMethod
Named method or configured method object.
**overrides : Any
Method fields overridden for this fit.
Returns
-------
TMLEResult or LongitudinalResult
Fitted result matching the study design.
Raises
------
CapabilityError
If the method is unavailable or the effect has no bound study.
MethodConfigurationError
If ``method`` does not satisfy the estimation-method contract.
See Also
--------
IdentifiedEffect.available_methods : Report supported methods before fitting.
TMLEMethod : Configure the default estimation method.
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")
... ),
... )
>>> effect = study.identify(ATE())
>>> result = effect.estimate(
... outcome_learner=LinearRegression(),
... treatment_learner=LogisticRegression(max_iter=1000),
... cross_fit=False,
... random_state=0,
... )
>>> sorted(result.estimates)
['ate']
"""
if self._study is None:
raise CapabilityError(
"this effect is not bound to analysis data; "
"construct a new CausalStudy to estimate it again",
reason=RefusalReason.INVALID_DECLARATION,
)
normalized = self._method(method, overrides)
refusal = self._unused_ratio_refusal(normalized)
if refusal is not None:
raise MethodConfigurationError.from_refusal(refusal)
if self.functional.longitudinal:
return self._estimate_longitudinal(normalized)
return self._estimate_point(normalized)
def _unused_ratio_refusal(self, method: EstimationMethod) -> Refusal | None:
"""Why ``ratio="classifier"`` would change nothing on this fit, or ``None``.
The setting selects how a policy's density ratio is estimated. A fit that
estimates no such ratio refuses it before any learner, rather than ignore it. A
point fit estimates one for ``policies=`` only, and not with a declared mechanism,
which fixes the ratio. A longitudinal fit estimates one at a continuous or vector
node only.
"""
models = getattr(method, "models", None)
if getattr(models, "ratio", "density") != "classifier":
return None
assert self._study is not None
data = self._study.data
if isinstance(data, LongitudinalData):
kinds = {data.node_kind(time) for time in range(1, data.n_times + 1)}
if kinds & {"continuous", "mixed"}:
return None
return Refusal(
RefusalReason.INVALID_DECLARATION,
"ratio='classifier' sets how the policy ratio of a continuous or vector node "
"is estimated, and every node of this design is categorical: its ratio is "
"the fitted probability of the assigned level.",
alternative="Drop ratio=.",
)
if self.functional.axis != "policy":
return Refusal(
RefusalReason.INVALID_DECLARATION,
"ratio='classifier' sets how a modified treatment policy's density ratio is "
f"estimated, and {type(self.estimand).__name__} declares no policy.",
alternative="Drop ratio=.",
)
if data.known_treatment is not None:
return KNOWN_POLICY_CLASSIFIER_REFUSAL
return None
def _method(
self, method: str | EstimationMethod, overrides: dict[str, Any]
) -> EstimationMethod:
if isinstance(method, str):
availability = {item.name: item for item in self.available_methods()}
record = availability.get(method)
if record is None:
raise CapabilityError(
f"unknown estimation method {method!r}; declarations: {list(availability)}",
reason=RefusalReason.INVALID_DECLARATION,
)
if not record.available:
raise CapabilityError.from_refusal(
_method_refusal(method, type(self.estimand).__name__, record.refusal)
)
constructors: dict[str, type[TMLEMethod]] = {
"tmle": TMLEMethod,
"collaborative_tmle": CollaborativeTMLEMethod,
"drtmle": DRTMLEMethod,
"sdr": SDRMethod,
"ipsi": IPSIMethod,
"gcomp": GComputationMethod,
"ipw": IPWMethod,
"parametric_gformula": ParametricGFormulaMethod,
}
return constructors[method]().with_overrides(**overrides)
if not isinstance(method, EstimationMethod):
raise MethodConfigurationError("method must be a named preset or an EstimationMethod")
declared = {item.name: item for item in self.available_methods()}.get(method.name)
if declared is None or not declared.available:
raise CapabilityError.from_refusal(
_method_refusal(
method.name,
type(self.estimand).__name__,
None if declared is None else declared.refusal,
)
)
return method.with_overrides(**overrides)
def _estimate_point(self, method: EstimationMethod) -> TMLEResult:
assert self._study is not None and isinstance(self._study.data, CausalData)
functional = self.functional
kwargs = method.estimator_kwargs()
kwargs.update({"estimands": (functional.target,), "reference": functional.reference})
if functional.axis == "regime":
kwargs["interventions"] = functional.interventions
elif functional.axis == "learned_rule":
(kwargs["learned_rule"],) = functional.interventions
elif functional.axis == "policy":
kwargs["policies"] = functional.interventions
elif functional.axis == "ipsi":
kwargs["incremental"] = functional.interventions
elif functional.axis == "msm":
kwargs["msm"] = functional.msm
if isinstance(method, CollaborativeTMLEMethod):
estimator: Any = CTMLE(**kwargs)
elif isinstance(method, DRTMLEMethod):
estimator = DRTMLE(**kwargs)
else:
estimator = TMLE(**kwargs)
result_set = estimator.fit(self._study.data)
raw = (
result_set[functional.intermediate]
if functional.intermediate is not None
else result_set.single()
)
raw = self._select_point_parameters(raw, method)
return replace(
raw,
identified_effect=self,
provenance=(
None
if raw.provenance is None
else replace(
raw.provenance,
protocol_fingerprint=(
None if self.protocol is None else self.protocol.fingerprint
),
)
),
method=method,
parameter_keys=self._point_parameter_keys(raw),
)
def _select_point_parameters(self, result: TMLEResult, method: EstimationMethod) -> TMLEResult:
if (
not isinstance(self.estimand, CounterfactualMean)
or self.estimand.treatment is None
or self.functional.target != "ey"
):
return result
data = result.data
if self.estimand.treatment not in data.treatment_levels:
raise DataError(
f"treatment {self.estimand.treatment!r} is not available; "
f"choose from {list(data.treatment_levels)}"
)
# Compose the alias from the arm label the estimator itself used, not from the value
# the caller typed. A numeric arm arrives as a float level, so ``treatment=1`` built
# ``ey[1]`` while the fit reported ``ey[1.0]``; nothing matched, and the narrowing
# returned an empty result rather than raising. Resolving through ``arm_label`` is
# the same "compose the known names forward" rule the parameter keys follow.
code = float(data.treatment_levels.index(self.estimand.treatment))
alias = parameter_name("ey", arm=data.arm_label(code))
estimates = {
name: value
for name, value in result.estimates.items()
if name == alias or name.startswith(f"{alias}[")
}
if not estimates:
raise CleverlyError(
f"selecting {self.estimand.treatment!r} kept none of the estimator's "
f"parameters {list(result.estimates)}; expected {alias!r}"
)
return replace(
result,
estimates=estimates,
simultaneous=_narrow_bands(result, estimates, method),
bootstrap=_narrow_bootstrap(result.bootstrap, estimates),
)
def _point_parameter_keys(self, result: TMLEResult) -> dict[str, ParameterKey]:
data = result.data
target = self.functional.target
reference_code = result.config.reference_arm
reference = data.arm_label(reference_code) if self.functional.axis == "arm" else None
base: dict[str, ParameterKey] = {}
if self.functional.axis == "arm":
if target == "ey":
for code in data.arm_codes:
value = data.arm_label(code)
alias = parameter_name("ey", arm=value)
base[alias] = ParameterKey(alias, target, value=value)
elif target == "ey1":
base[target] = ParameterKey(target, target, value=data.arm_label(1.0))
elif target == "ey0":
base[target] = ParameterKey(target, target, value=data.arm_label(0.0))
elif target in {"ate", "att", "atc", "rr", "or"}:
for code in data.arm_codes:
if code == reference_code:
continue
value = data.arm_label(code)
alias = arm_alias(
target,
arm=value,
versus=reference,
collapse=data.is_binary_treatment,
)
base[alias] = ParameterKey(alias, target, value=value, reference=reference)
elif target in {"par", "paf"}:
alias = arm_alias(target, arm=reference, collapse=data.is_binary_treatment)
base[alias] = ParameterKey(alias, target, value=reference)
elif target == "ey_obs":
base[target] = ParameterKey(target, target)
else:
state = {
"regime": result.nuisance.regimes,
# The fold-local rule is carried as a one-regime set named by the rule.
"learned_rule": result.nuisance.regimes,
"policy": result.nuisance.policies,
"rr_tilt": result.nuisance.policies,
"ipsi": result.nuisance.incremental,
"msm": result.nuisance.msm,
}[self.functional.axis]
if state is None:
raise CleverlyError(f"the {self.functional.axis} fit lost its structured state")
if self.functional.axis == "msm":
assert isinstance(state, MSMSet)
names = state.terms
else:
assert isinstance(state, (RegimeSet, PolicySet, IPSISet))
names = state.names
reference_name = names[int(state.reference)] if hasattr(state, "reference") else None
if target.startswith("ey_"):
for name in names:
alias = parameter_name(target, arm=name)
base[alias] = ParameterKey(alias, target, self.functional.axis, value=name)
elif target.startswith("ate_"):
for name in names:
if name == reference_name:
continue
alias = parameter_name(target, arm=name, versus=reference_name)
base[alias] = ParameterKey(
alias, target, self.functional.axis, value=name, reference=reference_name
)
else:
for term in names:
alias = parameter_name("msm", arm=term)
base[alias] = ParameterKey(alias, target, "msm", term=term)
keys = dict(base)
if data.has_strata:
for code in range(data.n_strata):
stratum = data.strata_levels[code]
label = data.stratum_label(code)
for alias, key in base.items():
conditional = stratum_alias(alias, label)
keys[conditional] = replace(key, alias=conditional, stratum=stratum)
selected = {name: keys[name] for name in result.estimates if name in keys}
if set(selected) != set(result.estimates):
raise CleverlyError(
"structured parameter keys disagree with the point estimator output: "
f"{list(selected)} != {list(result.estimates)}"
)
return selected
def _estimate_longitudinal(self, method: EstimationMethod) -> LongitudinalResult:
assert self._study is not None and isinstance(self._study.data, LongitudinalData)
functional = self.functional
engine: type[Any] = (
SDR
if isinstance(method, SDRMethod)
else IPSI
if isinstance(method, IPSIMethod)
else SequentialGComputation
if isinstance(method, GComputationMethod)
else SequentialIPW
if isinstance(method, IPWMethod)
else ParametricGFormula
if isinstance(method, ParametricGFormulaMethod)
else LTMLE
)
estimator = engine(
functional.interventions,
reference=functional.reference,
horizons=functional.horizons,
msm=functional.msm,
**method.estimator_kwargs(longitudinal=True),
)
raw = estimator.fit(self._study.data)
prefixes = {
"ey_regime": ("ey_regimen[", "risk_regimen[", "cif_regimen["),
"ate_regime": ("ate_regimen[", "risk_difference["),
"ey_ipsi": ("ey_regimen[", "risk_regimen[", "cif_regimen["),
"ate_ipsi": ("ate_regimen[", "risk_difference["),
"msm_regimen": ("msm_regimen[",),
}[functional.target]
estimates = {
name: value for name, value in raw.estimates.items() if name.startswith(prefixes)
}
index = None
if raw.parameter_index is not None:
index = {name: raw.parameter_index[name] for name in estimates}
raw = replace(
raw,
estimates=estimates,
parameter_index=index,
simultaneous=_narrow_bands(raw, estimates, method),
bootstrap=_narrow_bootstrap(raw.bootstrap, estimates),
)
return replace(
raw,
identified_effect=self,
provenance=replace(
raw.provenance,
protocol_fingerprint=(None if self.protocol is None else self.protocol.fingerprint),
),
method=method,
parameter_keys=self._longitudinal_parameter_keys(raw),
)
def _longitudinal_parameter_keys(self, result: LongitudinalResult) -> dict[str, ParameterKey]:
"""Structured keys read from the index the engine composed, not rebuilt beside it.
``LongitudinalResult.parameter_index`` already holds ``(label, cause, horizon)`` for
every reported parameter, recorded where the name was built. This used to rebuild that
grid as a cartesian product and zip it against the estimates positionally, so the
length guard could catch a missing cell but never a mis-*pairing*: swapping the cause
and horizon loops relabels every parameter of a competing-risks fit and raises nothing.
The reference was recovered by looking for ``"ate_"`` in the alias, which is the
display-name parsing the result contract forbids.
The index's label is the regimen for a mean and the engine's composed contrast label
for a contrast, so both spellings are composed forward from the declared regimens and
matched -- never split apart again.
"""
keys: dict[str, ParameterKey] = {}
if result.parameter_index:
reference = result.config.reference
regimens: dict[str, tuple[str, str | None]] = {}
for item in result.config.regimens:
regimens[item.label] = (item.label, None)
regimens[f"{item.label} vs {reference}"] = (item.label, reference)
for alias in result.estimates:
if alias not in result.parameter_index:
raise CleverlyError(
f"the longitudinal estimator reported {alias!r} with no entry in its "
"parameter index, so its regimen, cause, and horizon are unknown"
)
label, cause, horizon = result.parameter_index[alias]
if label not in regimens:
raise CleverlyError(
f"parameter {alias!r} is indexed by {label!r}, which is neither a "
f"declared regimen nor a contrast against {reference!r}"
)
regimen, against = regimens[label]
keys[alias] = ParameterKey(
alias,
self.functional.target,
"regimen",
value=regimen,
reference=against,
regimen=regimen,
cause=cause,
horizon=horizon,
)
elif result.msm is not None:
# A working-model fit has no parameter index -- its parameters are indexed by
# term -- so the names are composed forward from the declared terms and causes and
# matched. The previous rule took the term at ``position % len(terms)``, which
# assumed terms cycle fastest and dropped the cause entirely, so every parameter
# of a competing-risks projection claimed to be about no cause in particular.
composed = {
(f"msm_regimen[{term}]" if cause is None else f"msm_regimen[{term}, {cause}]"): (
term,
cause,
)
for cause in (result.config.causes or (None,))
for term in result.msm.terms
}
for name in result.estimates:
if name not in composed:
raise CleverlyError(
f"working-model parameter {name!r} is not composed from the declared "
f"terms {list(result.msm.terms)}, so its term is unknown"
)
term, cause = composed[name]
keys[name] = ParameterKey(name, "msm_regimen", "msm", term=term, cause=cause)
if set(keys) != set(result.estimates):
raise CleverlyError(
"structured parameter keys disagree with the longitudinal estimator output"
)
return keys