`statspai.tmle`¶

tmle ¶

Targeted Maximum Likelihood Estimation (TMLE) with Super Learner.

TMLE is a doubly robust, semiparametrically efficient estimator for causal effects that combines initial outcome regression with a targeted bias-correction step using the propensity score.

Components

TMLE : Full TMLE estimator for ATE/ATT with targeting step
SuperLearner : Ensemble learner for nuisance parameter estimation

References

van der Laan, M. J. & Rose, S. (2011). Targeted Learning: Causal Inference for Observational and Experimental Data. Springer Series in Statistics. [@vanderlaan2011targeted]

van der Laan, M. J., Polley, E. C., & Hubbard, A. E. (2007). Super Learner. Statistical Applications in Genetics and Molecular Biology, 6(1). [@vanderlaan2007super]

HALClassifier ¶

Bases: _BaseHAL

L1-penalised HAL logistic classifier (sklearn-compatible duck-typed API).

Fits an L1-penalised logistic regression on a main-effects HAL basis of per-feature step functions. Exposes .fit / .predict / .predict_proba so it can serve as the propensity learner in :func:sp.tmle and :func:sp.hal_tmle.

Parameters:

Name	Type	Description	Default
`C`	`float`	Inverse L1 penalty (larger = less shrinkage), as in scikit-learn.	`1.0`
`max_anchors_per_col`	`int`	Cap on step-function anchor points per feature.	`40`
`random_state`	`int`		`0`

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(200, 3))
>>> y = (X[:, 0] + rng.normal(size=200) > 0).astype(int)
>>> clf = sp.HALClassifier(max_anchors_per_col=10).fit(X, y)
>>> clf.predict_proba(X).shape
(200, 2)
>>> [int(c) for c in clf.classes_]
[0, 1]

HALRegressor ¶

Bases: _BaseHAL

L1-penalised HAL regressor (sklearn-compatible duck-typed API).

Fits an L1-penalised regression on a main-effects HAL basis of per-feature step functions. Exposes the .fit / .predict interface so it can be dropped into cross-fitting pipelines such as :func:sp.tmle and :func:sp.hal_tmle.

Parameters:

Name	Type	Description	Default
`lambda_`	`float`	L1 penalty. `None` selects it via cross-validation.	`None`
`max_anchors_per_col`	`int`	Cap on step-function anchor points per feature.	`40`
`cv`	`int`	Folds for the CV penalty search (used only when `lambda_` is None).	`5`
`random_state`	`int`		`0`

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(200, 3))
>>> y = X[:, 0] + 0.5 * X[:, 1] + rng.normal(size=200)
>>> reg = sp.HALRegressor(max_anchors_per_col=10).fit(X, y)
>>> reg.predict(X).shape
(200,)

LTMLEResult `dataclass` ¶

Bases: ResultProtocolMixin

Structured output of :func:ltmle.

Holds the treated/control marginal means, their ATE contrast, and the associated inference (se, ci, pvalue).

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 80
>>> l0 = rng.normal(size=n)
>>> a0 = rng.binomial(1, 1 / (1 + np.exp(-0.4 * l0)))
>>> l1 = 0.3 * l0 + 0.2 * a0 + rng.normal(size=n)
>>> a1 = rng.binomial(1, 1 / (1 + np.exp(-0.4 * l1)))
>>> y = 1.0 + 0.4 * a0 + 0.3 * a1 + 0.2 * l1 + rng.normal(scale=0.2, size=n)
>>> df = pd.DataFrame({"L0": l0, "A0": a0, "L1": l1, "A1": a1, "Y": y})
>>> res = sp.ltmle(
...     df, y="Y", treatments=["A0", "A1"],
...     covariates_time=[["L0"], ["L1"]],
... )
>>> isinstance(res, sp.LTMLEResult)
True
>>> float(res.ate)
0.71

LTMLEMSMResult `dataclass` ¶

Bases: ResultProtocolMixin

Working marginal structural model fitted by :func:ltmle_msm.

Attributes:

Name	Type	Description
`coefficients`	`DataFrame`	One row per working-model term: TMLE `estimate`, influence-curve `se`, `z`, `pvalue`, `ci_lower` / `ci_upper` and the IPTW estimate `iptw`.
`influence_curve`	`ndarray`	`n x p` influence curve of the TMLE coefficients.
`regime_means`	`DataFrame`	For each regime, the mean of the targeted `Q*_0(d)` (on the outcome's scale) and its weight.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n = 300
>>> L0 = rng.normal(size=n)
>>> A0 = rng.binomial(1, 0.5, size=n)
>>> L1 = 0.5 * L0 + rng.normal(size=n)
>>> A1 = rng.binomial(1, 0.5, size=n)
>>> Y = rng.binomial(1, 1 / (1 + np.exp(-(0.3 * (A0 + A1) + 0.2 * L1))))
>>> df = pd.DataFrame(dict(L0=L0, A0=A0, L1=L1, A1=A1, Y=Y))
>>> r = sp.ltmle_msm(df, 'Y', ['A0', 'A1'], [['L0'], ['L1']],
...                  [(0, 0), (1, 0), (0, 1), (1, 1)],
...                  {'cum_a': [0, 1, 1, 2]})
>>> isinstance(r, sp.LTMLEMSMResult)
True

LTMLESurvivalResult `dataclass` ¶

Bases: ResultProtocolMixin

Counterfactual survival curves and contrasts.

Produced by :func:ltmle_survival. Holds the treated/control counterfactual survival curves, the restricted-mean survival time (RMST) contrast with its standard error and confidence interval, and the terminal risk difference. Use :meth:to_frame for a tidy per-interval survival table.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n, K = 200, 3
>>> A = rng.integers(0, 2, size=(n, K))
>>> L = rng.normal(0, 1, size=(n, K))
>>> T = np.zeros((n, K), dtype=int)
>>> for k in range(K):
...     haz = 1 / (1 + np.exp(-(-1.5 - 0.5 * A[:, k] + 0.3 * L[:, k])))
...     T[:, k] = (rng.uniform(size=n) < haz).astype(int)
>>> df = pd.DataFrame({
...     **{f"T{k+1}": T[:, k] for k in range(K)},
...     **{f"A{k+1}": A[:, k] for k in range(K)},
...     **{f"L{k+1}": L[:, k] for k in range(K)},
... })
>>> res = sp.ltmle_survival(
...     df,
...     event_indicators=["T1", "T2", "T3"],
...     treatments=["A1", "A2", "A3"],
...     covariates_time=[["L1"], ["L2"], ["L3"]],
... )
>>> isinstance(res, sp.LTMLESurvivalResult)
True
>>> len(res.times)
3

SuperLearner ¶

Super Learner ensemble (van der Laan et al. 2007).

Parameters:

Name	Type	Description	Default
`library`	`list of sklearn estimators`	Candidate learners. If None, uses a default diverse library.	`None`
`n_folds`	`int`	Number of cross-validation folds.	`5`
`task`	`str`	'regression' or 'classification'.	`'regression'`
`random_state`	`int`		`42`

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> from sklearn.linear_model import LinearRegression, Ridge
>>> rng = np.random.default_rng(42)
>>> n = 400
>>> X = rng.normal(size=(n, 3))
>>> y = X @ np.array([1.0, -0.5, 0.3]) + rng.normal(scale=0.5, size=n)
>>> sl = sp.SuperLearner(
...     library=[LinearRegression(), Ridge(alpha=1.0)],
...     n_folds=3,
... ).fit(X, y)
>>> round(float(sl.weights_.sum()), 2)  # simplex weights
1.0
>>> sl.predict(X[:5]).shape
(5,)

References

vanderlaan2007super

get_params ¶

get_params(deep: bool = True) -> dict

Constructor parameters (scikit-learn get_params).

set_params ¶

set_params(**params: Any) -> SuperLearner

Set constructor parameters (scikit-learn set_params).

fit ¶

fit(X: Any, y: Any, sample_weight: Optional[Any] = None) -> SuperLearner

Fit the Super Learner.

Get cross-validated predictions from each base learner.
Find optimal weights via simplex-constrained least squares.
Refit all base learners on full data.

sample_weight (observation weights, e.g. survey weights) is passed to every base learner's fit and weights the simplex least-squares and the CV risks -- R SuperLearner's obsWeights. A learner whose fit does not take sample_weight raises rather than being fitted unweighted -- except in the default library, where such learners (the k-nearest neighbours learner) are left out with a warning and listed in dropped_learners_.

predict ¶

predict(X: Any) -> ndarray

Predict using the weighted ensemble.

For classification (task='classification') the returned values are the convex combination of base-learner probability predictions and are clipped to (1e-6, 1 - 1e-6) so that callers can take logit(.) without inf. For regression no clipping is applied.

Parameters:

Name	Type	Description	Default
`X`	`ndarray(n, p)`		required

Returns:

Type	Description
`ndarray(n)`

predict_proba ¶

predict_proba(X: Any) -> ndarray

Predict P(Y = 1) (for classification task).

Identical to :meth:predict under task='classification': a vector of shape (n,), not scikit-learn's (n, 2) array. The estimators that take a learner (sp.aipw, sp.ipw, sp.ltmle, sp.lmtp) accept both shapes.

summary ¶

summary() -> str

Print Super Learner summary.

TMLE ¶

Targeted Maximum Likelihood Estimation.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`		required
`y`	`str`		required
`treat`	`str`		required
`covariates`	`list of str`		required
`outcome_library`	`list of sklearn estimators`		`None`
`propensity_library`	`list of sklearn estimators`		`None`
`n_folds`	`int`		`5`
`estimand`	`str`		`'ATE'`
`alpha`	`float`		`0.05`
`propensity_bounds`	`tuple`		`(0.025, 0.975)`
`random_state`	`int`		`42`
`fold_indices`	`array-like of int`	Opt-in cross-validated TMLE partition (labels `0..K-1`, one per row of `data`); see :func:`statspai.tmle`. `None` (default) fits Q and g on the full sample, i.e. not cross-fitted.	`None`

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> x1 = rng.normal(size=n)
>>> x2 = rng.normal(size=n)
>>> ps = 1 / (1 + np.exp(-(0.5 * x1 - 0.3 * x2)))
>>> treat = rng.binomial(1, ps)
>>> y = 1.0 * treat + 0.8 * x1 - 0.5 * x2 + rng.normal(size=n)
>>> df = pd.DataFrame({'y': y, 'treat': treat, 'x1': x1, 'x2': x2})
>>> est = sp.TMLE(df, y='y', treat='treat', covariates=['x1', 'x2'],
...               n_folds=2, random_state=0)
>>> res = est.fit()
>>> bool(hasattr(res, 'estimate'))
True
>>> bool(res.se > 0)
True

fit ¶

fit() -> CausalResult

Run TMLE and return causal effect estimates.

ctmle ¶

ctmle(data: DataFrame, y: str, treat: str, covariates: List[str], estimand: str = 'ATE', outcome_library: 'Optional[List[BaseEstimator]]' = None, Q: 'Optional[np.ndarray]' = None, propensity_covariates: Optional[List[str]] = None, order: Optional[List[str]] = None, cv_folds: int = 5, fold_indices: 'Optional[Any]' = None, penalty: str = 'variance+bias', n_folds: int = 5, propensity_bounds: Tuple[float, float] = (0.025, 0.975), q_bound: float = 0.0005, alpha: float = 0.05, random_state: int = 42, cluster: Optional[str] = None, se_method: str = 'influence', n_boot: int = 200) -> CausalResult

Collaborative TMLE: the propensity model is chosen for the outcome fit.

Use it when a full propensity model produces extreme weights. A covariate enters the propensity only if targeting with it improves the fit of the outcome regression, so instruments and covariates the outcome model already accounts for tend to stay out, while a confounder the outcome model missed is brought in.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`		required
`y`	`str`	Outcome, binary or continuous (a continuous outcome is rescaled to the unit interval for the logistic fluctuation).	required
`treat`	`str`	Binary treatment (0/1).	required
`covariates`	`list of str`	Covariates of the outcome regression and, unless `propensity_covariates` says otherwise, the candidates for the propensity model.	required
`estimand`	`'ATE'`	Only the average treatment effect. For the effect on the treated use :func:`sp.tmle`; see Notes.	`'ATE'`
`outcome_library`	`list of sklearn estimators`	Super Learner library for the initial outcome regression `Q(A, W)`; the default library of :func:`sp.tmle` if omitted.	`None`
`Q`	`array of shape (n, 2)`	Initial outcome predictions `[Q(0, W), Q(1, W)]` on the outcome scale, replacing the Super Learner. They are then held fixed in the cross-validation, which is honest only if they were not fitted on these rows (or come from a parametric model with few parameters).	`None`
`propensity_covariates`	`list of str`	Candidates for the propensity model. Default: `covariates`.	`None`
`order`	`list of str`	A fixed order in which the candidates are offered to the propensity model, replacing the greedy search (the pre-ordered, "scalable" C-TMLE of Ju et al.). Each step then costs one propensity fit instead of one per remaining candidate, so the whole sequence is linear in the number of candidates where the greedy search is quadratic. Candidates not listed are never offered. Put first what is most likely to confound.	`None`
`cv_folds`	`int`	Folds of the cross-validation that selects the step.	`5`
`fold_indices`	`array - like`	One fold label per row of `data`, replacing the random assignment (and `cv_folds`).	`None`
`penalty`	`('variance+bias', 'variance', 'search', 'none')`	What is added to the residual sum of squares when candidates are compared. The variance of the estimator's influence function is what keeps a covariate that only makes the weights extreme out of the propensity model. `'search'`: the variance is added in the greedy search and the cross-validation compares residual sums of squares alone. This is what R `ctmle::ctmleDiscrete` does by default. `'variance'`: the variance is added in the cross-validation too. `'variance+bias'`: the cross-validation also adds `n` times the squared mean of the held-out influence function, the estimator's squared bias, as in Gruber and van der Laan (2010). `'none'`: residual sums of squares throughout.	`'variance+bias'`
`n_folds`	`int`	Internal folds of the Super Learner.	`5`
`propensity_bounds`	`tuple`	Truncation of every fitted propensity.	`(0.025, 0.975)`
`q_bound`	`float`	Truncation of the initial outcome predictions on the unit scale.	`5e-4`
`alpha`	`float`		`0.05`
`random_state`	`int`	Seeds the fold assignment and the Super Learner.	`42`
`cluster`	`str`	Cluster column. Folds are formed from whole clusters and the standard error sums the influence function within clusters.	`None`
`se_method`	`('influence', 'bootstrap')`	`'bootstrap'` reruns the whole procedure, selection included, on `n_boot` resamples of rows (of clusters with `cluster=`); the standard error is the standard deviation of the estimates and the interval their percentile interval. Prefer it for inference: see Notes. Not with `Q`.	`'influence'`
`n_boot`	`int`		`200`

Returns:

Type Description

CausalResult

detail has one row per step of the sequence: the covariate added, the estimate at that step, the empirical and the cross-validated criterion with its three parts, and which step was selected. model_info holds selected_covariates (in the order they entered), candidate_order (the full greedy order), step, n_restarts and the per-row influence_function. After se_method='bootstrap' it also has selection_frequency, the share of resamples in which each candidate entered the propensity model.

Notes

What is selected and how. The propensity models are logistic regressions on subsets of the candidates, nested along a greedy order. Step 0 is an intercept-only propensity. At each step the covariate that gives the lowest criterion is added; when none lowers it, the current targeted fit becomes the initial fit and the search goes on. The step is then chosen by cross-validation: the sequence is rebuilt on each training fold and its steps are scored on the held-out rows.

Inference. The influence-function standard error is the standard deviation of the efficient influence function at the selected targeted fit and propensity model, over sqrt(n). It ignores the selection, and when the selected propensity model is not the true one (which is the point of the method) that function is not the estimator's influence function. In the simulations behind this function it was 20% to 40% too small and its 95% interval covered 61% to 88% of the time; R ctmle behaves the same way. se_method='bootstrap' is the remedy.

When it helps and when it does not. With an outcome model that is about right and a covariate that drives treatment but not the outcome, the root mean squared error was 27% below that of :func:sp.tmle with the full propensity model, and the bootstrap interval covered 94% of the time. With an outcome model that omits a strong confounder the confounder is brought into the propensity first, as it should be, but the error was no smaller than that of :func:sp.tmle. It is a tool for weak overlap, not a default.

Why only the ATE. A collaborative ATT was implemented and tested. With an additive outcome model and effects that vary with a confounder, the intercept-only propensity already gives the lowest criterion, nothing rewards a step that would repair the missing heterogeneity, and the estimate was biased by 0.17 to 0.20 where :func:sp.tmle with the full propensity was unbiased. R ctmle estimates the ATE only as well.

What it does not protect against. A confounder that neither the outcome model nor the candidate list contains.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> from sklearn.linear_model import LinearRegression
>>> rng = np.random.default_rng(0)
>>> n = 600
>>> x1, z = rng.normal(size=n), rng.normal(size=n)
>>> a = rng.binomial(1, 1 / (1 + np.exp(-(0.5 * x1 + 2.0 * z))))
>>> yv = a + x1 + rng.normal(size=n)          # z is an instrument
>>> df = pd.DataFrame({'y': yv, 'a': a, 'x1': x1, 'z': z})
>>> res = sp.ctmle(df, y='y', treat='a', covariates=['x1', 'z'],
...                outcome_library=[LinearRegression()])
>>> 'z' in res.model_info['selected_covariates']
False

References

[@vanderlaan2010collaborative], [@gruber2010application], [@ju2019scalable]

statspai.tmle¶

tmle ¶

HALClassifier ¶

HALRegressor ¶

LTMLEResult dataclass ¶

LTMLEMSMResult dataclass ¶

LTMLESurvivalResult dataclass ¶

SuperLearner ¶

get_params ¶

set_params ¶

fit ¶

predict ¶

predict_proba ¶

summary ¶

TMLE ¶

fit ¶

ctmle ¶

`statspai.tmle`¶

LTMLEResult `dataclass` ¶

LTMLEMSMResult `dataclass` ¶

LTMLESurvivalResult `dataclass` ¶