r"""Simulation studies: does the estimator work, and does its interval cover?
Everything else in this package validates a *fit*. This validates the *estimator*:
repeat the whole procedure on fresh data from a process whose truth is known, and
measure what comes out.
The three numbers that matter, in order:
**Coverage.** The fraction of replications whose confidence interval contains the
truth. A nominal 95% interval that covers 80% of the time is not conservative or
approximate -- it is wrong, and every p-value derived from it is wrong too. This is
the single most informative check available, and it is the one most often skipped.
**Bias, scaled by root-n.** Absolute bias can shrink while remaining large relative to
the standard-error scale. :attr:`StudyResult.root_n_bias` therefore reports
:math:`\sqrt{n}\,\mathrm{bias}`. Bounded values support a negligible first-order bias
claim; they do not establish efficiency. Growth suggests that a residual bias term is
not vanishing fast enough.
**Standard-error accuracy.** The ratio of the mean estimated standard error to the
actual standard deviation of the estimates across replications. Near 1 means the
reported uncertainty is honest; below 1 means it is optimistic, which is how a
coverage failure usually arises.
This harness is public API, not just test scaffolding: running it on a DGP resembling
your own data is the most direct way to find out whether a particular
estimator configuration is trustworthy for your problem before you rely on it.
"""
from __future__ import annotations
import os
import pickle
from collections import Counter
from collections.abc import Callable, Iterable, Sequence
from dataclasses import dataclass
from pathlib import Path
from typing import TYPE_CHECKING, Any, Protocol, cast, runtime_checkable
import numpy as np
from .._inference_status import (
STATUS_RECORDS,
InferenceStatus,
is_labelled,
precedent_status,
status_record,
supplies_inference,
)
from .._typing import FloatArray
from ..estimators.base import TMLEResultSet
from ..estimators.direct_effect import check_level
from ..inference.influence import spread_name
from ..utils.parallel import map_parallel
from ..utils.text import format_table
from ._failures import failure_call, recordable
from ._seeds import sample_seed_streams, validate_sample_seeds
if TYPE_CHECKING: # pragma: no cover - typing only
from ..datasets.synthetic import DGP
from ..estimators.tmle import TMLE
from ..longitudinal import LTMLE
__all__ = [
"CoverageStudy",
"EstimandSummary",
"ReplicationFailure",
"ReplicationRecord",
"StudyResult",
"summarize_replications",
]
[docs]
@dataclass(frozen=True)
class ReplicationRecord:
"""One estimand from one successful repeated-sampling draw.
Parameters
----------
replicate : int
Index of the draw.
seed : int
Seed that produced it.
estimand : str
Alias this record is for.
truth : float
The estimand's true value under the process.
estimate : float
What the estimator reported.
std_error : float
The standard error it reported.
covered : bool
Whether the interval contained ``truth``.
rejected : bool
Whether the null was rejected at ``alpha``.
inference_estimate : float
The estimate on the inference scale.
alpha : float
Significance level the interval was built at.
inference : str
Which quantity ``std_error``, ``covered`` and ``rejected`` measured. One of
:data:`~cleverly.inference.influence.InferenceStatus`, which the
:doc:`inference reference </technical-reference/inference>` lists.
``"influence_curve"`` says the estimator's own reported inference. Any other
status says the estimator supplies none, and the study measured its retained
plug-in diagnostic instead.
support : {"supported", "experimental"}
The support label of the replicate's estimate.
"""
replicate: int
seed: int
estimand: str
truth: float
estimate: float
std_error: float
covered: bool
rejected: bool
inference_estimate: float
alpha: float
inference: InferenceStatus = "influence_curve"
support: str = "experimental"
[docs]
@dataclass(frozen=True)
class ReplicationFailure:
"""A replication whose fit failed numerically, kept so the study can be diagnosed.
Parameters
----------
replicate : int
Index of the draw that failed.
seed : int
Seed that produced it, so the failure can be reproduced.
error_type : str
Exception class raised.
message : str
Its message.
call : str
``module:function`` of the frame that raised it.
"""
replicate: int
seed: int
error_type: str
message: str
call: str = ""
ReplicationOutcome = tuple[ReplicationRecord, ...] | ReplicationFailure
[docs]
@runtime_checkable
class ReplicationCheckpoint(Protocol): # numpydoc ignore=PR01
"""Storage that lets an interrupted :class:`CoverageStudy` resume.
:meth:`CoverageStudy.run` calls ``load`` for every replication before it fits any, and
fits only those that return ``None``. It calls ``save`` in the worker process as soon as
a replication finishes, so a stopped study loses only the fits in flight. When
``n_jobs`` is not 1, ``save`` therefore runs in another process: an implementation must be
picklable and must store outside the process, as :class:`DirectoryCheckpoint` does.
See Also
--------
DirectoryCheckpoint : The implementation that keeps one file per replication.
"""
[docs]
def load(self, replicate: int) -> ReplicationOutcome | None:
"""Return the stored outcome of ``replicate``, or ``None`` when there is none.
Parameters
----------
replicate : int
Index of the replication.
Returns
-------
tuple of ReplicationRecord, ReplicationFailure or None
The stored outcome, or ``None``.
"""
...
[docs]
def save(self, replicate: int, outcome: ReplicationOutcome) -> None:
"""Store the outcome of ``replicate``.
Parameters
----------
replicate : int
Index of the replication.
outcome : tuple of ReplicationRecord or ReplicationFailure
What the replication returned.
"""
...
[docs]
@dataclass(frozen=True)
class DirectoryCheckpoint:
"""Keep each finished replication of a :class:`CoverageStudy` as a file in a directory.
One file per replication, written under a temporary name and then renamed, so a study
stopped mid-write never reads a partial file. The directory belongs to one study. The
load compares each stored seed with the study's seed and refuses a mismatch. It compares
nothing else, so a study with the same seed and another ``n``, data-generating process or
estimator reads the stored results. Give each configuration its own directory.
Parameters
----------
directory : str or os.PathLike
Where the files go. It is created when missing.
See Also
--------
ReplicationCheckpoint : The protocol this class implements.
Examples
--------
>>> import tempfile
>>> from sklearn.linear_model import LinearRegression, LogisticRegression
>>> from cleverly.datasets import nonlinear_dgp
>>> from cleverly.estimators import TMLE
>>> from cleverly.validation import CoverageStudy, DirectoryCheckpoint
>>> def study(directory):
... return CoverageStudy(
... dgp=nonlinear_dgp(),
... estimator=lambda: TMLE(
... outcome_learner=LinearRegression(),
... treatment_learner=LogisticRegression(),
... cross_fit=False,
... ),
... n=200,
... n_replicates=3,
... seed=0,
... checkpoint=DirectoryCheckpoint(directory),
... )
>>> with tempfile.TemporaryDirectory() as directory:
... first = study(directory).run()
... resumed = study(directory).run()
>>> resumed.replications == first.replications
True
"""
directory: str | os.PathLike[str]
def _path(self, replicate: int) -> Path:
return Path(self.directory) / f"replication-{replicate:09d}.pkl"
[docs]
def load(self, replicate: int) -> ReplicationOutcome | None:
"""Return the stored outcome of ``replicate``, or ``None`` when there is none.
Parameters
----------
replicate : int
Index of the replication.
Returns
-------
tuple of ReplicationRecord, ReplicationFailure or None
The stored outcome, or ``None``.
"""
path = self._path(replicate)
if not path.exists():
return None
with path.open("rb") as handle:
return cast("ReplicationOutcome", pickle.load(handle))
[docs]
def save(self, replicate: int, outcome: ReplicationOutcome) -> None:
"""Store the outcome of ``replicate``.
Parameters
----------
replicate : int
Index of the replication.
outcome : tuple of ReplicationRecord or ReplicationFailure
What the replication returned.
"""
path = self._path(replicate)
path.parent.mkdir(parents=True, exist_ok=True)
temporary = path.with_suffix(f".{os.getpid()}.tmp")
with temporary.open("wb") as handle:
pickle.dump(outcome, handle, protocol=pickle.HIGHEST_PROTOCOL)
os.replace(temporary, path)
@dataclass(frozen=True)
class EstimandSummary:
"""Simulation results for one estimand."""
estimand: str
truth: float
n: int
n_replicates: int
estimates: FloatArray
std_errors: FloatArray
covered: FloatArray
rejected: FloatArray
#: The point estimates on the **inference scale** -- the scale ``std_errors`` is on, which
#: for a ratio estimand is the log scale and for everything else is ``estimates`` itself.
#: ``None`` means the two coincide, which is the case for every difference and level.
#: Only :attr:`se_ratio` reads it; see there for why it has to exist.
inference_estimates: FloatArray | None = None
#: What ``std_errors``, ``covered`` and ``rejected`` measured: the estimator's own
#: inference, or its retained plug-in diagnostic when it supplies none.
#: :func:`summarize_replications` reads it off the records. When the replicates took
#: more than one status, it is the one :func:`precedent_status` gives.
inference: InferenceStatus = "influence_curve"
#: The statuses the replicates took, each with its replicate count, when they took
#: more than one; empty when every replicate took :attr:`inference`. A fit whose
#: status depends on the draw, such as a clustered fit whose cluster count straddles
#: the few-cluster floor, gives a mix.
status_counts: tuple[tuple[str, int], ...] = ()
#: ``"experimental"`` when any replicate's estimate was, else ``"supported"``.
support: str = "experimental"
@property
def mixed_statuses(self) -> str:
"""The replicates' statuses and counts as one phrase, or ``""`` when uniform."""
return ", ".join(f"{status} in {count}" for status, count in self.status_counts)
@property
def mean_estimate(self) -> float:
return float(np.mean(self.estimates))
@property
def inference_scale_estimates(self) -> FloatArray:
""":attr:`estimates` on the scale the reported standard error is on."""
if self.inference_estimates is None:
return self.estimates
return self.inference_estimates
@property
def bias(self) -> float:
return float(np.mean(self.estimates) - self.truth)
@property
def root_n_bias(self) -> float:
r"""``sqrt(n) * bias``: bounded as ``n`` grows for a root-n consistent estimator."""
return float(np.sqrt(self.n) * self.bias)
@property
def monte_carlo_se(self) -> float:
"""Standard deviation of the estimates across replications: the *actual* variability."""
return float(np.std(self.estimates, ddof=1))
@property
def bias_se(self) -> float:
"""Standard error of the estimated bias, for judging whether it is real."""
return float(self.monte_carlo_se / np.sqrt(self.n_replicates))
@property
def mean_std_error(self) -> float:
"""Mean *reported* standard error."""
return float(np.mean(self.std_errors))
@property
def se_ratio(self) -> float:
"""Reported over actual standard error; 1.0 is honest, below 1 is optimistic.
Both halves are taken on the **inference scale**, which is the scale
:attr:`~cleverly.inference.ParameterEstimate.std_error` is defined on -- the log scale
for a ratio, and the reporting scale for everything else. For a ratio the two come
apart badly: ``psi`` is the odds ratio while ``std_error`` is ``SE(log OR)``, and by
the delta method ``sd(OR) ~ psi * sd(log OR)``, so dividing one by the other returns
roughly ``1 / psi`` and says nothing whatever about calibration. A well-behaved
``or`` of ``0.42`` came back as ``2.82``, which read as an interval three times too
wide and was only the reciprocal of the truth.
:attr:`monte_carlo_se` deliberately stays on the reporting scale, because that is
where :attr:`bias` and :attr:`rmse` are meaningful; the two differ only for a ratio.
"""
spread = float(np.std(self.inference_scale_estimates, ddof=1))
if spread <= 0:
return float("nan")
return float(self.mean_std_error / spread)
@property
def coverage(self) -> float:
"""Fraction of confidence intervals containing the truth."""
return float(np.mean(self.covered))
@property
def coverage_se(self) -> float:
"""Standard error of the coverage estimate -- needed to judge a near-miss."""
p = self.coverage
return float(np.sqrt(p * (1.0 - p) / self.n_replicates))
@property
def rejection_rate(self) -> float:
"""Fraction of replications rejecting the null.
Power when the truth is non-zero; the type I error rate when it is zero.
"""
return float(np.mean(self.rejected))
@property
def rmse(self) -> float:
return float(np.sqrt(np.mean((self.estimates - self.truth) ** 2)))
def to_dict(self) -> dict[str, Any]:
# The convention ``ParameterEstimate.to_dict`` keeps: an ordinary summary carries
# the columns it always carried, and a summary of a diagnostic carries
# ``inference`` and names its mean spread through ``spread_name``.
row: dict[str, Any] = {"estimand": self.estimand}
if is_labelled(self.inference):
row["inference"] = self.inference
if self.status_counts:
row["mixed_statuses"] = self.mixed_statuses
if self.support == "experimental":
row["support"] = self.support
return {
**row,
"truth": self.truth,
"n": self.n,
"n_replicates": self.n_replicates,
"mean_estimate": self.mean_estimate,
"bias": self.bias,
"bias_se": self.bias_se,
"root_n_bias": self.root_n_bias,
"rmse": self.rmse,
"monte_carlo_se": self.monte_carlo_se,
spread_name("mean_std_error", self.inference): self.mean_std_error,
"se_ratio": self.se_ratio,
"coverage": self.coverage,
"coverage_se": self.coverage_se,
"rejection_rate": self.rejection_rate,
}
def summarize_replications(
records: Sequence[ReplicationRecord], *, estimand: str, n: int
) -> EstimandSummary:
"""Build the canonical descriptive summary for one estimand's records.
Every record measures the plug-in numbers of its replicate: the interval and standard
error an inferential fit reports are the same numbers a diagnostic fit retains. So
records that took different statuses summarize together, under the status
:func:`~cleverly._inference_status.precedent_status` gives, and the summary names the
statuses it read in :attr:`EstimandSummary.status_counts`.
"""
selected = tuple(record for record in records if record.estimand == estimand)
if not selected:
raise ValueError(f"no successful replication records for estimand {estimand!r}")
estimates = np.asarray([record.estimate for record in selected], dtype=float)
inference = np.asarray([record.inference_estimate for record in selected], dtype=float)
counts: Counter[str] = Counter(record.inference for record in selected)
status = precedent_status(counts)
# The table order, then "influence_curve", which has no record and comes last.
order = [*STATUS_RECORDS, "influence_curve"]
return EstimandSummary(
estimand=estimand,
truth=float(np.mean([record.truth for record in selected])),
n=n,
n_replicates=len(selected),
estimates=estimates,
std_errors=np.asarray([record.std_error for record in selected], dtype=float),
covered=np.asarray([record.covered for record in selected], dtype=float),
rejected=np.asarray([record.rejected for record in selected], dtype=float),
inference_estimates=None if np.array_equal(inference, estimates) else inference,
inference=status,
support=(
"experimental"
if any(record.support == "experimental" for record in selected)
else "supported"
),
status_counts=(
tuple((name, counts[name]) for name in order if name in counts)
if len(counts) > 1
else ()
),
)
def _diagnostic_noun(statuses: Iterable[str]) -> str:
"""Use a specific noun only when every diagnostic has the same status."""
diagnostic = {status for status in statuses if not supplies_inference(status)}
if len(diagnostic) != 1:
return "plug-in diagnostic"
return status_record(next(iter(diagnostic))).diagnostic_noun
[docs]
@dataclass(frozen=True)
class StudyResult:
"""The full output of a :class:`CoverageStudy`.
Parameters
----------
summaries : dict of str to EstimandSummary
Bias, coverage and rejection rate per estimand.
replications : tuple of ReplicationRecord
One record per estimand per successful draw.
failures : tuple of ReplicationFailure
One record per draw that raised.
n : int
Sample size per replication.
n_replicates : int
Replications requested.
alpha : float
Significance level coverage was measured at.
label : str
Name of the study, for reports.
"""
summaries: dict[str, EstimandSummary]
replications: tuple[ReplicationRecord, ...]
failures: tuple[ReplicationFailure, ...]
n: int
n_replicates: int
alpha: float
label: str
@property
def n_failed(self) -> int:
"""Return the number of failed replications."""
return len(self.failures)
def __getitem__(self, estimand: str) -> EstimandSummary:
return self.summaries[estimand]
[docs]
def to_frame(self, backend: str | None = None) -> Any:
"""Return tabular output in the input dataframe backend.
Parameters
----------
backend : {"pandas", "polars", "pyarrow"} or None, default=None
Dataframe backend to return. ``None`` uses pandas when installed, then the first
available backend.
Returns
-------
dataframe
One row per estimand, with bias and coverage.
"""
from ..utils.frames import frame_from_dict
rows = [summary.to_dict() for summary in self.summaries.values()]
# The union of the keys, in first-seen order, as ``TMLEResult.to_frame`` builds
# it: a summary of a diagnostic carries columns an ordinary one does not.
keys = list(dict.fromkeys(key for row in rows for key in row))
payload = {key: [row.get(key) for row in rows] for key in keys}
return frame_from_dict(payload, backend=backend)
@property
def _diagnostic(self) -> bool:
"""Whether any summary measured a plug-in diagnostic, not inference."""
return any(not supplies_inference(summary.inference) for summary in self.summaries.values())
@property
def _diagnostic_noun(self) -> str:
"""What the diagnostic columns measured across all summaries and replicates."""
return _diagnostic_noun(
status
for summary in self.summaries.values()
for status, _ in (summary.status_counts or ((summary.inference, summary.n_replicates),))
)
[docs]
def summary(self) -> str:
"""Return a printable summary.
Returns
-------
str
A printable table, one line per estimand.
"""
level = 1.0 - self.alpha
lines = [
f"Simulation study: {self.label}",
"-" * (18 + len(self.label)),
f"n = {self.n}; replications = {self.n_replicates}"
+ (f" ({self.n_failed} failed)" if self.n_failed else ""),
# Printed beside the budget rather than in a column, because it describes
# every row: it says whether "mean se" and the coverage column measured the
# estimator's own inference or its retained diagnostic.
*(
[
f"measuring a {self._diagnostic_noun}: one or more reported "
"estimates supply no interval"
+ (
" on some replicates"
if any(
any(status == "influence_curve" for status, _ in summary.status_counts)
for summary in self.summaries.values()
)
else ""
)
+ ", so the columns below describe the spread of "
"the curve it reports and not a confidence interval"
]
if self._diagnostic
else []
),
*(
f"{summary.estimand}: the replicates took more than one status "
f"({summary.mixed_statuses}), and every column reads them as the "
f"{_diagnostic_noun(status for status, _ in summary.status_counts)}"
if any(not supplies_inference(status) for status, _ in summary.status_counts)
else f"{summary.estimand}: the replicates took more than one status "
f"({summary.mixed_statuses}), and every column reads the interval each "
"replicate reported"
for summary in self.summaries.values()
if summary.status_counts
),
*(
f"{summary.estimand}: experimental, because at least one replicate's "
"estimate has no registered study behind its configuration"
for summary in self.summaries.values()
if summary.support == "experimental"
),
"",
format_table(
[
"estimand",
"truth",
"mean",
"bias",
"sqrt(n)*bias",
"mc se",
"mean se",
"se ratio",
f"{level:.0%} cover",
"reject",
],
[
[
summary.estimand,
f"{summary.truth:.4f}",
f"{summary.mean_estimate:.4f}",
f"{summary.bias:+.4f}",
f"{summary.root_n_bias:+.3f}",
f"{summary.monte_carlo_se:.4f}",
f"{summary.mean_std_error:.4f}",
f"{summary.se_ratio:.3f}",
f"{summary.coverage:.3f}",
f"{summary.rejection_rate:.3f}",
]
for summary in self.summaries.values()
],
),
"",
self.verdict(),
]
return "\n".join(lines)
[docs]
def verdict(self) -> str:
"""Reading of the study, naming the specific failure mode where there is one.
Returns
-------
str
A reading of the study, naming the failure mode where there is one.
"""
notes: list[str] = []
target = 1.0 - self.alpha
for summary in self.summaries.values():
shortfall = target - summary.coverage
if shortfall > 3.0 * summary.coverage_se and shortfall > 0.02:
notes.append(
f"{summary.estimand}: coverage {summary.coverage:.3f} is below the nominal "
f"{target:.0%} by more than Monte Carlo error"
+ (
f"; the {spread_name('reported standard error', summary.inference)} is "
f"{1 / summary.se_ratio:.2f}x too small"
if summary.se_ratio < 0.95
else "; the standard error looks right, so this is bias"
)
)
if abs(summary.bias) > 3.0 * summary.bias_se and abs(summary.bias) > 0.02 * max(
1.0, abs(summary.truth)
):
notes.append(
f"{summary.estimand}: bias {summary.bias:+.4f} is "
f"{abs(summary.bias) / summary.bias_se:.1f} Monte Carlo standard errors from "
"zero"
)
if not notes and self._diagnostic:
# A plug-in interval that covers is not a confidence interval the estimator
# reports, so the verdict does not certify one. Bias is a property of the
# point estimate, which this path does report.
return (
"VERDICT: bias is consistent with a correctly working estimator. The "
f"coverage column measures a {self._diagnostic_noun}, so it "
"certifies no confidence interval for this estimator."
)
if not notes:
return "VERDICT: coverage and bias are consistent with a correctly working estimator."
return "VERDICT:\n" + "\n".join(f" - {note}" for note in notes)
def __repr__(self) -> str: # pragma: no cover - debugging aid
return self.summary()
[docs]
class CoverageStudy:
"""Repeatedly fit an estimator on fresh data and measure bias and coverage.
Parameters
----------
dgp : DGP or callable
A :class:`~cleverly.datasets.DGP`, or a callable
``(n, seed) -> (frame, truth)`` following the convention of the generators in
:mod:`cleverly.datasets`.
estimator : callable
A zero-argument factory returning a fresh estimator per replication. A factory
rather than an instance, so replications cannot share fitted state.
n : int
Sample size per replication.
n_replicates : int
How many replications to run.
estimands : sequence of str or None
Which estimands to summarise; defaults to whatever the first fit reports.
fit_kwargs : mapping or None
Passed to ``fit``. Column names, ``delta=``, ``id=`` and so on.
seed : int or None
Root of the distinct per-replication sample seeds.
sample_seeds : sequence of int or None
Explicit distinct uint32 sample seeds, one per replication. Use these to coordinate
separate studies. When supplied, these replace the sample stream from ``seed``.
n_jobs : int
Number of joblib workers across replications.
truth_key : {"population", "sample"}
``"population"`` (default) compares against the population estimand, fixed
across replications; ``"sample"`` compares against each replication's realised
sample estimand, which removes one source of variability but changes what
coverage means.
intermediate_value : object or None
The level of the intermediate variable to study, for a controlled direct effect.
Required when ``fit_kwargs`` contains ``intermediate=``, because such a fit
returns one result per level and each level is a *different parameter*: the
study picks that level out of the result set and asks the process for the truth
at the same level, so the two cannot silently disagree. See
:mod:`cleverly.estimators.direct_effect`.
label : str
Name of the study, used in its reports.
checkpoint : ReplicationCheckpoint or None
Storage for finished replications. With one, :meth:`run` reuses every stored
outcome and fits only the rest, so a stopped study resumes where it stopped.
Notes
-----
A replication whose fit raises a numerical or fitting error is recorded as a
:class:`ReplicationFailure` and the study continues. These are recorded:
:class:`~cleverly.exceptions.DataError`, :class:`~cleverly.exceptions.LongitudinalError`,
``numpy.linalg.LinAlgError``, ``FloatingPointError``, ``ZeroDivisionError``,
``OverflowError``, and a plain ``ValueError`` or ``RuntimeError`` that ``numpy``, ``scipy``,
``scikit-learn`` or ``statsmodels`` raised. Every other exception stops the study,
including :class:`~cleverly.exceptions.CapabilityError`,
:class:`~cleverly.exceptions.MethodConfigurationError` and programming errors such as
``TypeError``. An error while drawing a sample always stops the study.
Examples
--------
>>> from sklearn.linear_model import LinearRegression, LogisticRegression
>>> from cleverly.estimators import TMLE
>>> from cleverly.datasets import nonlinear_dgp
>>> from cleverly.validation import CoverageStudy
>>> study = CoverageStudy(
... dgp=nonlinear_dgp(),
... estimator=lambda: TMLE(
... outcome_learner=LinearRegression(),
... treatment_learner=LogisticRegression(max_iter=1000),
... ),
... n=500,
... n_replicates=50,
... seed=0,
... )
>>> study.n_replicates
50
"""
def __init__(
self,
dgp: DGP | Callable[..., tuple[Any, dict[str, float]]],
estimator: Callable[[], TMLE | LTMLE],
*,
n: int = 1000,
n_replicates: int = 100,
estimands: Sequence[str] | None = None,
fit_kwargs: dict[str, Any] | None = None,
seed: int | None = None,
sample_seeds: Sequence[int] | None = None,
n_jobs: int = 1,
truth_key: str = "population",
intermediate_value: float | None = None,
label: str | None = None,
checkpoint: ReplicationCheckpoint | None = None,
) -> None:
if n_replicates < 2:
raise ValueError(f"n_replicates must be at least 2; got {n_replicates}")
if truth_key not in ("population", "sample"):
raise ValueError(f"truth_key must be 'population' or 'sample'; got {truth_key!r}")
if intermediate_value is not None:
intermediate_value = check_level(intermediate_value)
# Checked here rather than in the replication loop, which swallows every
# exception to keep one bad draw from killing a study: a mismatch between the
# fit and the level asked for is a configuration error, and reporting it as
# "every replication failed" would bury it.
targets_a_level = (fit_kwargs or {}).get("intermediate") is not None
if targets_a_level and intermediate_value is None:
raise ValueError(
"fit_kwargs names an intermediate variable, so each fit returns one "
"result per level and each level is a different parameter. Pass "
"intermediate_value=0.0 or 1.0 to say which one this study measures."
)
if intermediate_value is not None and not targets_a_level:
raise ValueError(
f"intermediate_value={intermediate_value} was given but fit_kwargs does "
"not name an intermediate variable, so the fit has no level to select."
)
self.dgp = dgp
self.estimator = estimator
self.n = n
self.n_replicates = n_replicates
self.estimands = estimands
self.fit_kwargs = fit_kwargs or {"outcome": "Y", "treatment": "A"}
self.seed = seed
self.sample_seeds = (
None if sample_seeds is None else validate_sample_seeds(sample_seeds, n_replicates)
)
self.n_jobs = n_jobs
self.truth_key = truth_key
self.intermediate_value = intermediate_value
self.label = str(label or getattr(dgp, "name", getattr(dgp, "__name__", "study")))
self.checkpoint = checkpoint
def _draw(self, seed: int) -> tuple[Any, dict[str, float]]:
from ..datasets.synthetic import DGP as DGPClass
if isinstance(self.dgp, DGPClass):
# Passing the level through is what keeps the comparison honest: the truth a
# process reports for a controlled direct effect is the effect *at a level*,
# and DGP.sample silently defaults it to 0 when it is not told which one.
return self.dgp.sample(self.n, seed=seed, intermediate_value=self.intermediate_value)
if self.intermediate_value is None:
return self.dgp(self.n, seed=seed)
return self.dgp(self.n, seed=seed, intermediate_value=self.intermediate_value)
def _select(self, result: Any) -> Any:
"""Pick the single result a replication is summarising.
A fit with ``intermediate=`` returns one result per level, and the levels are
different parameters rather than two views of one, so the study has to be told
which it is measuring coverage for rather than guessing. That is the set's own
key lookup: ``intermediate_value`` is ``None`` for an ordinary fit, which is
exactly the key such a fit uses, so a mismatch in either direction surfaces as a
``KeyError`` naming the levels that are available.
An estimator that does not return a *set* -- :class:`~cleverly.LTMLE` returns one
:class:`~cleverly.LongitudinalResult` -- is already the single result, and is
passed through. Keying into it would ask for the parameter named ``None`` and
get a ``KeyError`` that the replication loop swallows, so the whole study would
fail with "every replication failed; check the estimator configuration" while the
estimator was configured correctly.
"""
if not isinstance(result, TMLEResultSet):
if self.intermediate_value is not None:
raise TypeError(
f"intermediate_value={self.intermediate_value!r} was given, but the "
f"estimator returned a {type(result).__name__} rather than a result "
"set with one entry per level"
)
return result
return result[self.intermediate_value]
[docs]
def run_replication(
self,
replicate_index: int,
seed: int,
) -> tuple[ReplicationRecord, ...] | ReplicationFailure:
"""Fit one declared replication from an explicit sample seed.
A fresh estimator is created for every call. A recordable fitting failure (see the
class notes) returns an explicit record; any other exception propagates. This method
does not aggregate a sampling distribution or change the study budget.
Parameters
----------
replicate_index : int
Original replication index within the declared budget.
seed : int
Explicit uint32 seed supplied to the data generator.
Returns
-------
tuple of ReplicationRecord or ReplicationFailure
One record per estimand on success, or the fit's explicit failure record.
"""
import warnings
if (
isinstance(replicate_index, bool)
or not isinstance(replicate_index, int)
or not 0 <= replicate_index < self.n_replicates
):
raise ValueError("replicate_index must be an integer within the declared budget")
seed = validate_sample_seeds((seed,), 1)[0]
# The draw stays outside the handler: an error there is a defect of the process,
# not a fit that broke down on one sample.
frame, truth = self._draw(int(seed))
try:
with warnings.catch_warnings():
# Individual replications routinely trip positivity warnings; the
# aggregate coverage is the diagnostic here, not the per-fit warnings.
warnings.simplefilter("ignore")
fitted = self.estimator().fit(frame, **self.fit_kwargs)
result = self._select(fitted)
names = self.estimands or tuple(result.estimates)
out: list[ReplicationRecord] = []
for name in names:
estimate = result[name]
prefix = "" if self.truth_key == "population" else "sample_"
reference = truth[f"{prefix}{name}"]
# This class is the instrument that *measures* the calibration of
# whatever a fit reports, so it reads the retained diagnostic when the
# estimator supplies no inference. Refusing here would make it
# impossible to measure the diagnostic the package deliberately keeps,
# and, worse, the `except Exception` below would convert the refusal
# into a `ReplicationFailure`: every replication of every selector cell
# would be dropped silently and the cell would publish on a shrunken
# budget. `inference` on the record says which one was measured.
# No branch on the status: the plug-in accessors read the private
# bodies ``ci`` and ``std_error`` read, and ``_plugin_pvalue`` is the
# body ``pvalue`` returns, so an inferential estimate gives the numbers
# it always gave and a diagnostic one gives the refused properties'
# arithmetic rather than a second rule that could disagree with it.
low, high = estimate.plugin_interval
error = estimate.plugin_std_error
probability = estimate._plugin_pvalue()
# The fifth entry is the estimate on the scale `std_error` is on, which is
# the log scale for a ratio -- taken off `log_psi`, the same field `ci`
# builds the interval from, rather than re-derived here. `se_ratio` is the
# only consumer, and it is the only summary that has to compare the two.
out.append(
ReplicationRecord(
replicate=replicate_index,
seed=int(seed),
estimand=name,
truth=float(reference),
estimate=float(estimate.psi),
std_error=float(error),
covered=bool(low <= reference <= high),
rejected=bool(probability < estimate.alpha),
inference_estimate=estimate.inference_value,
alpha=float(estimate.alpha),
inference=estimate.inference,
support=estimate.support.label,
)
)
return tuple(out)
except Exception as error:
if not recordable(error):
raise
return ReplicationFailure(
replicate=replicate_index,
seed=int(seed),
error_type=type(error).__name__,
message=str(error),
call=failure_call(error),
)
def _checkpointed(self, replicate_index: int, seed: int) -> ReplicationOutcome:
"""One replication, stored in the checkpoint as soon as it finishes."""
outcome = self.run_replication(replicate_index, seed)
if self.checkpoint is not None:
self.checkpoint.save(replicate_index, outcome)
return outcome
[docs]
def run(self) -> StudyResult:
"""Execute the study.
Returns
-------
StudyResult
Per-replication records, per-estimand summaries, and the failures.
"""
if self.sample_seeds is None:
root = (
self.seed if self.seed is not None else cast(int, np.random.SeedSequence().entropy)
)
seeds = sample_seed_streams({root: self.n_replicates})[root]
else:
seeds = validate_sample_seeds(self.sample_seeds, self.n_replicates)
stored: dict[int, ReplicationOutcome] = {}
if self.checkpoint is not None:
for index, seed in enumerate(seeds):
outcome = self.checkpoint.load(index)
if outcome is None:
continue
if isinstance(outcome, ReplicationFailure):
kept: int | None = outcome.seed
else:
kept = outcome[0].seed if outcome else None
if kept is not None and kept != int(seed):
raise RuntimeError(
f"the checkpoint holds replication {index} under seed {kept}, but "
f"this study draws it from seed {seed}"
)
stored[index] = outcome
missing = [(index, seed) for index, seed in enumerate(seeds) if index not in stored]
fitted = dict(
zip(
(index for index, _ in missing),
map_parallel(self._checkpointed, missing, n_jobs=self.n_jobs),
strict=True,
)
)
outcomes = [
stored[index] if index in stored else fitted[index] for index in range(len(seeds))
]
failures = tuple(outcome for outcome in outcomes if isinstance(outcome, ReplicationFailure))
records = tuple(
record
for outcome in outcomes
if not isinstance(outcome, ReplicationFailure)
for record in outcome
)
if not records:
detail = ""
if failures:
first = failures[0]
detail = f"; first failure was {first.error_type}: {first.message}"
raise RuntimeError(
"every replication failed; check the estimator configuration" + detail
)
alphas = {record.alpha for record in records}
if len(alphas) != 1:
raise RuntimeError(f"replications returned inconsistent alpha levels: {sorted(alphas)}")
alpha = alphas.pop()
names = tuple(dict.fromkeys(record.estimand for record in records))
summaries = {
name: summarize_replications(records, estimand=name, n=self.n) for name in names
}
successful_draws = len({record.replicate for record in records})
return StudyResult(
summaries=summaries,
replications=records,
failures=failures,
n=self.n,
n_replicates=successful_draws,
alpha=alpha,
label=self.label,
)