Longitudinal treatment#
Declare temporal roles#
LongitudinalTreatment aligns treatment nodes, time-varying histories, censoring indicators, and
outcome processes. At treatment node t, a dynamic rule sees only history available by t.
For the reported influence-curve inference, the rule must be fixed before the fit and rowwise. A
callable may vectorize over the history frame, but it must not estimate a threshold from the
sample or otherwise make one row’s assignment depend on other rows. The rule needs a declaration,
which Declare a dynamic rule shows.
from sklearn.linear_model import LinearRegression, LogisticRegression
from cleverly import CausalStudy, LongitudinalTreatment, RegimeContrast
from cleverly.datasets import make_longitudinal
frame, truth = make_longitudinal(n=2_000, seed=11)
study = CausalStudy(
frame,
design=LongitudinalTreatment(
outcome="Y",
treatment=("A1", "A2"),
baseline=("W1", "W2"),
time_varying=((), ("L2",)),
censoring=("C1", "C2"),
),
)
result = study.estimate(
RegimeContrast({"always": 1, "never": 0}, reference="always"),
outcome_learner=LinearRegression(),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
n_folds=3,
learner_folds=3,
random_state=0,
)
The resolved regimen matrix is shared by nuisance fitting, follower masks, targeting, and report keys. A regimen cannot look ahead or change interpretation between stages. An empirically learned policy needs inference for the policy-learning step and is outside this estimator’s contract.
With n_folds > 1, each outer training fold runs an untargeted backward regression sequence. One
pooled fluctuation per node then targets the out-of-fold predictions over every follower. The fit
stores the realized assignment in result.folds.
Cross-fitting the recursion
gives the steps and the published result they follow.
Each node’s fluctuation solves its score over every follower, as on a single-fold fit.
res.diagnostics.score_equations() reports one solver row per node.
The cross-fitted end-of-study study
measured standard-error ratios from 0.9797 to 1.0130 over 1,600 replications at n=2000. Those
results apply only to that study’s law and estimator settings. They do not establish calibrated
variance for other laws, weights, clusters, survival outcomes, or sample sizes.
Read a long panel#
Customer data usually arrives as one row per customer and decision time. Each row holds the
covariates known before the decision. It also holds the decision and the observation indicator.
LongPanel and LongitudinalData.from_long_panel reshape such a frame into the trajectory
container that LongitudinalTreatment builds from a wide frame. The reshape ends in
LongitudinalData.from_frame, so every check of the wide layout runs on the result.
argument |
the column holds |
wide name |
|---|---|---|
|
the customer identifier. One value is one trajectory |
the row of the container, kept in |
|
the decision time. Values must sort |
the node index |
|
the decision at that row. A list of columns is one vector node |
|
|
unit-level covariates, measured before the first decision. Name at least one column |
the original name |
|
covariates measured at that row, before its decision |
|
|
|
|
|
an end-of-study outcome, one value per customer, on any of its rows |
the original name |
|
a survival event indicator on each row. A mapping of cause to column declares competing risks |
|
|
every decision time, in time order. Omit it when the times are numbers, dates or durations, and the reader sorts them |
the node order |
|
an account or store identifier, constant within a customer |
|
|
a unit-level probability weight |
|
|
the known mechanism, as |
|
Declare exactly one of outcome= and event=. One customer trajectory is one observation for the
inference model. Repeated rows of one customer add no independent observations. Declare
cluster= when customers depend on each other, such as the customers of one account.
A generated name is "column@t", where t is the node index and not the time label. The labels
stay in data.panel, a PanelLayout, which data.panel.decision_times lists. A rule at node t
reads "column@s" for s <= t. The summary of a fit on an adapted panel adds a long panel:
line. The line names the unit column, the decision times, and any column left out at a node.
from_long_panel takes rows in any order. The container follows the sorted unit identifiers.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression, LogisticRegression
from cleverly import CausalStudy, LongPanel, RegimeContrast
rng = np.random.default_rng(5)
rows, bought = [], {}
for customer in range(400):
age, eng, offers = rng.normal(), rng.normal(), 0
for week in (10, 20, 30):
offer = rng.random() < 1 / (1 + np.exp(-0.5 * eng))
stays = week == 30 or rng.random() < 0.9
rows.append((customer, week, age, eng, "offer" if offer else "none", float(stays)))
if not stays:
break
offers += offer
eng = 0.6 * eng + 0.5 * offer + rng.normal(scale=0.7)
else:
p = 1 / (1 + np.exp(-(0.5 * eng + 0.4 * offers - 0.5)))
bought[customer] = float(rng.random() < p)
panel = pd.DataFrame(rows, columns=["customer", "week", "age", "eng", "offer", "active"])
panel["bought"] = panel["customer"].map(bought)
long_design = LongPanel(
unit="customer",
time="week",
treatment="offer",
baseline=["age"],
time_varying=["eng"],
censoring="active",
outcome="bought",
)
long_study = CausalStudy(panel, design=long_design)
print(long_study.data.history_names(2), long_study.data.panel.decision_times)
The panel holds 400 customers over weeks 10, 20 and 30. A customer who leaves has active equal to
0 on its last row, and no outcome. history_names(2) returns ('age', 'eng@1', 'eng@2'). The
design hands the frame to CausalStudy, and from_long_panel takes the same arguments if you
want the container first. The sequential offers tutorial
reads a long panel with a three-level offer and declared offer probabilities.
Put each measurement on the right row#
A time_varying column is measured before that row’s decision. Code cannot read the order of a
measurement from the data, so the declaration states it. A measurement taken after the decision,
such as the response to the offer, belongs to the next row. Move it with one line of pandas:
after = panel.assign(response=panel["eng"]) # recorded after each row's decision
after["response"] = after.groupby("customer")["response"].shift()
The first row of each customer then holds a missing value. A leading run of missing values is not
an error. The column starts at its first recorded node, and data.panel.omitted lists each column
and node left out. A column that is missing at a later node, after it was first recorded, is
refused. Leaving it out there would remove a confounder from that node’s adjustment set.
A column that holds the previous treatment is not a leak. Among the units that follow a plan, it
equals the arm that the plan assigned. Under a policy, the earlier observed arm is what
policy_frame already exposes to a density.
Write one rule for every node#
each_node(rule, n_times) returns one rule per node. Node t calls rule(history, t), so the
rule reads f"eng@{t}" and nothing later. Put the result in a DynamicRegimen, as
Declare a dynamic rule states.
from cleverly.longitudinal import DynamicRegimen, each_node
engaged = DynamicRegimen(
"engaged",
each_node(lambda history, t: np.where(history[f"eng@{t}"] > 0, "offer", "none"), 3),
rule_kind="known",
)
long_result = long_study.estimate(
RegimeContrast({"never": ("none",) * 3, "engaged": engaged}, reference="never"),
outcome_learner=LogisticRegression(max_iter=1000),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
censoring_learner=LogisticRegression(max_iter=1000),
n_folds=2,
random_state=0,
)
print(long_result.summary())
This plan offers at node t when the customer’s engagement at node t is positive. The fit
raises DataError for a rule that names a column from a later node. It raises the same error when
the column is absent at the node, for example because the node screened it out as constant. It
never reads an earlier node instead.
test |
what it checks |
|---|---|
|
|
|
a rule raises when its column is absent at its node |
|
|
Go back to a long panel#
data.to_long_panel() returns one row per customer and node at which the customer is still in the
study. Rows are sorted by unit and time, in the backend that the data came from.
back = long_study.data.to_long_panel()
print(back.head(3))
The result differs from the input in three ways.
difference |
reason |
|---|---|
the weights have mean one |
the container normalizes the weights |
a covariate screened out at a node is missing at that node |
the node holds no column for it |
the known probability columns are not returned |
the container keeps them as mechanisms, not as covariates |
to_long_panel raises DataError on a container built from a wide frame.
What the reader refuses#
Each refusal raises DataError before any learner. A message about one unit names the first such unit, and the decision time when one applies.
the panel holds |
why the reader refuses it |
|---|---|
time labels that are not numbers, dates or durations, with no |
the sort order of a label such as |
a missing unit or time, a duplicated |
a long panel holds one row per unit and decision time. |
a customer whose first row is after the first decision, or a gap between two rows |
late entry is a different observation scheme. Intermittent observation is not censoring |
rows after a customer left, or rows that stop with no censoring 0 and no event |
a customer who is censored or has the event leaves for good. Declare |
an event value on a row whose censoring is 0 |
the order at a node is \(A_t, C_t, Y_t\). Record an event as censoring 1 and event 1. Record a loss after an event-free interval as censoring 0 and a missing event |
a |
a later value would let a rule at an earlier node read the future |
an |
an end-of-study outcome holds one value per unit. A per-row event indicator is |
|
a long panel takes one outcome column. Declare a repeated outcome on a wide frame |
one column in two roles, or a covariate that equals the outcome or the row’s event |
a rule could read the outcome under another name |
a time-varying column with no value at a node after its first recorded node |
the column would leave that node’s adjustment set |
|
the package fits no intercept-only nuisance design. The refusal has code |
a treatment that never changes within a customer |
that is one baseline decision held over every node. Build the wide frame and declare |
a column named like a generated name, such as |
the name collides with a generated column. Rename it |
The wide layout refuses the same event on a censored row. LongitudinalData.from_frame used to
drop it without a message, and it now raises DataError.
Long horizons#
The registered evidence for LTMLE and SDR covers two decision nodes. The
roadmap plans laws with
three to five. Every \(L_t\) block accumulates into the history design, so a long panel adds columns
at each node and needs more customers. times= cannot coarsen a grid. To model fewer decisions,
aggregate the frame to the coarser decision times before you read it.
Declare a dynamic rule#
Write each plan that holds a rule as a DynamicRegimen, and declare rule_kind="known". The
declaration states that every callable node is a fixed rowwise function of the history.
from cleverly.longitudinal import DynamicRegimen
adaptive = DynamicRegimen(
"adaptive",
(1, lambda history: (history["L2"] > 0).astype(int)),
rule_kind="known",
)
rule_result = study.estimate(
RegimeContrast({"never": 0, "adaptive": adaptive}, reference="never"),
outcome_learner=LinearRegression(),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
n_folds=3,
learner_folds=3,
random_state=0,
)
This plan treats at the first node. At the second node it treats each unit whose L2 is positive.
DynamicRegimen raises CapabilityError for rule_kind=None and for "estimated". The fit
refuses a callable written inline in regimens=, such as {"adaptive": (1, rule)}. An inline
callable carries no declaration. A plan of labels alone, such as {"never": 0}, needs no
declaration. Code cannot inspect a closure, so the fit accepts a false "known" declaration.
Declare a stochastic policy#
A node can draw its arm from a known distribution. Write the node as a Stochastic node with
density_kind="known", inside a DynamicRegimen. The density function receives the history and
the earlier treatments, and it returns one column per treatment level, in sorted level order.
import numpy as np
from cleverly.interventions import Stochastic
def half(history):
return np.full((len(history), 2), 0.5)
def by_marker(history):
treat = np.where(history["L2"] > 0, 0.75, 0.25)
return np.column_stack([1.0 - treat, treat])
mostly = DynamicRegimen(
"mostly",
(
Stochastic(half, "half", density_kind="known"),
Stochastic(by_marker, "marker", density_kind="known"),
),
)
policy_result = study.estimate(
RegimeContrast({"never": 0, "mostly": mostly}, reference="never"),
outcome_learner=LinearRegression(),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
n_folds=3,
learner_folds=3,
random_state=0,
)
This plan treats half of the units at the first node. At the second node it treats three quarters
of the units whose L2 is positive and a quarter of the others. The parameter is the mean outcome
if every unit drew its arms this way.
A rule node cannot read an earlier treatment, because the history it receives holds none. To
continue the arm that an earlier policy node drew, write a one-hot Stochastic node that reads
it. The fit resolves a one-hot density as the rule it equals.
def continue_drawn(history):
drawn = (history["A1"] == 1).to_numpy(dtype=float)
return np.column_stack([1.0 - drawn, drawn])
continued = DynamicRegimen(
"continued",
(
Stochastic(half, "half", density_kind="known"),
Stochastic(continue_drawn, "continue", density_kind="known"),
),
)
Positivity is required wherever the policy puts mass. A policy that can draw a level that no unit
on the plan received at that node raises LongitudinalError before any learner.
Known stochastic policies
gives the estimator, its conditions, and its evidence. Code cannot inspect a closure, so the fit
accepts a false "known" declaration.
Survival outcomes#
An outcome sequence declares one absorbing event process and makes horizon part of the estimand.
from cleverly import RegimeMean
from cleverly.datasets import make_longitudinal_survival
frame, truth = make_longitudinal_survival(n=2_000, seed=12)
study = CausalStudy(
frame,
design=LongitudinalTreatment(
outcome=("Y1", "Y2"),
treatment=("A1", "A2"),
baseline=("W1", "W2"),
time_varying=((), ("L2",)),
censoring=("C1", "C2"),
),
)
risks = study.estimate(
RegimeMean({"always": 1, "never": 0}, horizons=(1, 2)),
outcome_learner=LinearRegression(),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
)
Each node fit uses the relevant uncensored, event-free risk set. The result key preserves regimen and horizon.
Read ratios and the restricted mean survival time (RMST) off the fitted curve. They refit nothing.
risk_ratio = risks.ratio("risk_regimen[always @ t=2]", "risk_regimen[never @ t=2]")
survival_ratio = risks.ratio(
"risk_regimen[always @ t=2]", "risk_regimen[never @ t=2]", view="survival"
)
rmst = risks.rmst("always", 3, versus="never")
print(risk_ratio.ci, survival_ratio.psi, rmst.psi, rmst.ci)
call |
what it reports |
|---|---|
|
the risk ratio of two levels at one horizon, with its interval on the log scale |
|
the ratio of the two survival probabilities |
|
the odds ratio of the two levels |
|
the RMST up to node |
|
the difference of two RMSTs |
|
the restricted mean time lost to one cause |
rmst needs the risk at every horizon from 1 to tau - 1. For nodes spaced \(\Delta\) apart,
multiply by \(\Delta\) to get calendar time.
Pass n_bootstrap= to refit the whole estimator on resamples. Each estimate then carries a
percentile interval beside its influence-curve interval.
bootstrapped = study.estimate(
RegimeMean({"always": 1, "never": 0}, horizons=(1, 2)),
outcome_learner=LinearRegression(),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
n_bootstrap=200,
random_state=0,
)
print(bootstrapped.summary())
The bootstrap contract
states what each replicate resamples and refits. The percentile interval carries the support label
of the fit. A registered support row decides the label. The bootstrap kind of the row covers an
in-sample end-of-study or survival fit with static regimens. Any other fit prints
[experimental].
The registered study measures correctly specified cell-mean nuisances on finite
binary laws, at \(n = 1000\). A fit with fewer than 20 clusters prints bootstrap sd and a
percentile range, a diagnostic. ratio and rmst carry the bootstrap of their inputs.
Bootstrap for sparse regimens#
This feature is experimental. The X39 long-horizon studies have not run, so every item in this
section prints experimental. The label makes no coverage claim. The roadmap row
X39 holds the pending
studies.
Use the refit bootstrap for LTMLE or SDR when few units follow a regimen, or when the
influence curve has heavy tails. The Wald interval reads its variance from those few rows. The
interval then under-covers, because the variance estimate is noisy and the sampling distribution is
skewed.
A Wald fit warns with SparseInfluenceWarning when this happens. A fit with a percentile
bootstrap warns too, because the percentile interval keeps the skew. A fit with BCa does not warn.
The bootstrap is opt in. A fit never turns it on.
option |
what it does |
label today |
|---|---|---|
|
refits the whole estimator on |
|
|
takes the quantiles of the |
as above |
|
reads bias-corrected and accelerated limits from the draws. It needs |
|
|
adds a studentized bootstrap sup-t band to the fit, for two or more estimates |
|
|
reports how few rows carry each influence curve. Call it directly: |
a diagnostic, with no label |
Pass these options to LTMLE or SDR. The CausalStudy.estimate route does not accept
bootstrap_interval.
from sklearn.linear_model import LinearRegression, LogisticRegression
from cleverly.datasets import make_longitudinal
from cleverly.longitudinal import LTMLE
frame, _ = make_longitudinal(n=300, seed=0)
result = LTMLE(
{"always": 1, "never": 0},
outcome_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
censoring_learner=LogisticRegression(max_iter=1000),
n_folds=1,
n_bootstrap=20,
bootstrap_interval="bca",
simultaneous=True,
random_state=0,
).fit(
frame,
outcome="Y",
treatment=["A1", "A2"],
baseline=["W1", "W2"],
time_varying=[[], ["L2"]],
censoring=["C1", "C2"],
)
print(result.summary())
print(result.bootstrap_band.critical_value)
The example uses 20 resamples to run fast. Use at least 199 in an analysis. With B resamples a
BCa limit is an order statistic of the draws, so a small B gives a coarse limit.
Percentile or BCa. The percentile interval takes the 2.5% and 97.5% quantiles of the draws. It
assumes the draws are symmetric about the estimate. The BCa interval shifts both limits by a bias
correction, from the share of draws below the estimate, and by an acceleration, from the skew of
the influence curve. Choose BCa when diagnostics.influence() reports a large skewness. For a
symmetric curve the two intervals are close.
The sup-t band. Pass simultaneous=True with n_bootstrap to get result.bootstrap_band.
Each replicate records its own standard error. The band takes the largest studentized deviation
across the estimates in each replicate. Its critical value is a high order statistic of those
maxima.
The band is \(\hat\psi_j \pm c^* \widehat{SE}_j\) for each estimate. The multiplier band in
result.simultaneous does not change. The summary lists both.
The influence report. result.diagnostics.influence() returns one row for each estimate.
column |
meaning |
|---|---|
|
the effective number of rows that carry the curve, \(n / \hat\kappa\). With |
|
the third and fourth standardized moments of the curve. A normal curve has kurtosis 3 |
|
the relative standard deviation of the variance estimate |
|
the Wald coverage shortfall that the second-order Edgeworth term predicts. A positive value is a loss and a negative value is a predicted excess. The flag does not read this column |
|
the BCa bias correction and acceleration. They are |
The warning fires when variance_rel_sd is above 0.125, which is close to a Kish count of 64. The
warning lists the relative SD of each flagged row, with its Kish count and Edgeworth loss. A
flagged row can show a negative loss, because the skew term can predict an excess while the noisy
variance estimate still triggers the flag. The
X39 work set that threshold after it read the pilot runs, so it is not a pre-registered rule.
A cross-fitted fit, a hazard-model fit and an isotonic fit warn too, but the message does not
advise the bootstrap. For the cross-fitted fit the bootstrap over-covers, and the roadmap owner
X20-bootstrap in the roadmap holds that result.
Cost. Each resample refits every nuisance. A fit of 64,000 rows and five decisions takes about
0.2 seconds for each resample on the fast path, a development measurement that no committed
artifact records yet. The X39 studies plan B = 199, so one such fit costs about 40 seconds of
bootstrap time.
The fast path fits the distinct rows with the resample counts as weights. It
applies to an in-sample, unclustered, unweighted fit with declared treatment and censoring
mechanisms, categorical treatment, and learners that declare frequency_exact. Any other fit
refits each materialized resample, which costs much more. Set n_jobs to run resamples in
parallel.
Evidence. The code ships as experimental. The long-horizon Python studies over three to five
decisions are pending, as a follow-up. The R comparators (ltmle and lmtp) are deferred to
X59. No coverage result for
BCa, the sup-t band or the warning threshold exists in this repository yet.
Repeated outcomes#
A repeated outcome is measured at every node and can return to zero. Declare it with
Repeated([...]) as the outcome. The fit reports the mean of the outcome at every horizon.
from cleverly import RegimeContrast, Repeated
from cleverly.datasets import make_longitudinal_repeated
frame, truth = make_longitudinal_repeated(n=1_500, seed=7)
study = CausalStudy(
frame,
design=LongitudinalTreatment(
outcome=Repeated(["Y1", "Y2", "Y3"]),
treatment=("A1", "A2", "A3"),
baseline=("W",),
time_varying=((), ("L2",), ("L3",)),
censoring=("C1", "C2", "C3"),
),
)
means = study.estimate(
RegimeContrast({"never": 0, "always": 1}, reference="never"),
outcome_learner=LogisticRegression(max_iter=1000),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
)
print(means.curve(scale="mean")[["estimand", "psi", "ci_lower", "ci_upper"]])
item |
value |
|---|---|
result keys |
|
|
one row for each reported key and horizon. It is the only scale |
survival view |
not available. A repeated outcome has no absorbing event |
missing |
refused. Encode the loss as |
|
refused before any learner. |
truth holds the exact mean at each horizon, under the names the fit reports. In the generator, W
takes three values and every later node is binary. The fit can therefore check its curve against
truth.
Repeated outcomes gives the
estimator, its conditions, and its evidence. The
repeated-outcome study
checks its intervals against exact truths. Its comparison with R ltmle 1.3-0 is pending.
The repeated-outcome IPSI study
does the same for IPSI. It has no comparator.
Competing risks#
A mapping from cause labels to absorbing outcome sequences declares competing risks. Cause, horizon, and regimen remain separate key fields. The engine uses a cause-specific recursion while keeping all prior competing events out of later risk sets.
A baseline treatment and a time-to-event outcome#
Use TimeToEvent when each unit has one follow-up time, one event code and one treatment given at
baseline. Code 0 means censored. The design bins the times onto a grid of visit times and holds
the treatment over every visit.
from cleverly import TimeToEvent
from cleverly.datasets import make_point_survival
frame, truth = make_point_survival(n=1_000, seed=4, n_times=4)
design = TimeToEvent(
time="time", event="event", treatment="A", baseline=("W1", "W2"), grid=(1, 2, 3, 4)
)
curve = CausalStudy(frame, design=design).estimate(
RegimeContrast({"arm1": 1, "arm0": 0}, reference="arm0", horizons=(2, 4)),
outcome_learner=LinearRegression(),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
cross_fit=False,
)
print(curve["ate_regimen[arm1 vs arm0 @ t=4]"].ci, truth["ate_regimen[arm1 vs arm0 @ t=4]"])
Follow these rules for the input.
rule |
reason |
|---|---|
Declare |
a grid read off the sample would change with the sample |
Record a censoring time as a grid time, or as a time after the last grid time |
a censoring between two grid times cannot be ordered against that interval’s events, so the fit refuses it |
Give |
the parameter name keeps the node index, and |
Expect \(K(K+1)/2\) regressions per regimen and cause for \(K\) grid times |
a default grid on day-level times can hold thousands of nodes |
Expect a risk of zero with no interval at a grid time before any follower of an arm had the event |
the regression of a node with no event is zero, its maximum-likelihood hazard. That node adds no variance to the curve, so the risk and its contrasts report no interval, and the fit warns. |
An event after the last grid time leaves the unit event-free at the end of the grid. causes=
maps each nonzero code to a cause label. Without it, codes 0 and 1 declare one event and any
other codes declare competing causes. Declare continuous_treatment=True for a continuous dose,
and give each regimen as a modified treatment policy such as Shift(0.5, cap=None).
On the wide layout, declare the same design with one column name: treatment="A" on
LongitudinalTreatment. The node count then comes from the outcome columns, else from
censoring=, else from time_varying=. The wide layout also takes a covariate measured after
baseline, and an end-of-study outcome after several censoring nodes. A regimen is one arm, a known
rule, a known stochastic policy, or a modified treatment policy, all at baseline. A plan that
changes the treatment after baseline is refused, because the design holds one decision.
One baseline treatment held over the nodes
gives the contract and its evidence.
A survival curve that never decreases#
The default survival fit targets each horizon in its own backward pass, so its curve can decrease by a small amount. Two options give a non-decreasing curve.
option |
what it does |
what it needs |
|---|---|---|
|
fits the cause-specific hazard at every node, and targets one set of hazards for every horizon together |
one categorical treatment at baseline, held over every node, and no covariate after baseline |
|
projects each fitted curve and its band limits onto a non-decreasing curve |
any survival or competing-risk fit |
Pass the hazard learner as outcome_learner. A hazard fit does not use pseudo_learner.
from sklearn.linear_model import LinearRegression, LogisticRegression
from cleverly import CausalStudy, RegimeMean, TimeToEvent
from cleverly.datasets import make_point_survival
frame, truth = make_point_survival(n=1_000, seed=4, n_times=4)
study = CausalStudy(
frame,
design=TimeToEvent(
time="time", event="event", treatment="A", baseline=("W1", "W2"), grid=(1, 2, 3, 4)
),
)
hazard_curve = study.estimate(
RegimeMean({"arm1": 1, "arm0": 0}),
survival_model="hazard",
targeting="one_step",
outcome_learner=LogisticRegression(max_iter=1000),
treatment_learner=LogisticRegression(max_iter=1000),
censoring_learner=LogisticRegression(max_iter=1000),
cross_fit=False,
)
print(hazard_curve.curve()[["estimand", "psi", "ci_lower", "ci_upper"]])
The risks of each arm do not decrease over the four grid times. Call summary() to see the walk: the
parameter count, the step count, the largest relative score, and whether the walk converged.
targeting="iterative", the default, targets the hazards at each horizon on its own. That curve is
not guaranteed to be monotone.
To keep the iterated-mean fit and project its curve, pass isotonic=True.
projected = study.estimate(
RegimeMean({"arm1": 1, "arm0": 0}),
isotonic=True,
outcome_learner=LinearRegression(),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(max_iter=1000),
cross_fit=False,
)
print(projected.isotonic_moves)
isotonic_moves lists each level the projection moved. On this sample it is empty, because the
fitted curve already does not decrease. The projection keeps every influence curve and standard
error. A pointwise interval at a moved level is valid only where the true curve strictly
increases. The isotonic projection
states each claim and its condition.
The hazard model refuses a treatment at several nodes, a covariate after baseline, a continuous treatment, and a modified treatment policy. Each message names the iterated-mean route. Hazard-based survival targeting gives the contract, both algorithms, and the evidence. No registered study measures the coverage of a hazard fit or a projected fit yet.
Fit the sequentially doubly robust estimator#
Pass method="sdr" to estimate() to fit the sequentially doubly robust estimator. It reports
the same parameters as the default method. It stays consistent when, at each node, either the
mechanism or the regression of the node’s transform is consistent.
from sklearn.linear_model import LinearRegression, LogisticRegression
from cleverly import CausalStudy, LongitudinalTreatment, RegimeContrast
from cleverly.datasets import make_longitudinal
frame, _ = make_longitudinal(n=300, seed=0)
study = CausalStudy(
frame,
design=LongitudinalTreatment(
outcome="Y",
treatment=("A1", "A2"),
baseline=("W1", "W2"),
time_varying=((), ("L2",)),
censoring=("C1", "C2"),
),
)
result = study.identify(RegimeContrast({"always": 1, "never": 0}, reference="never")).estimate(
method="sdr",
outcome_learner=LogisticRegression(),
pseudo_learner=LinearRegression(),
treatment_learner=LogisticRegression(),
censoring_learner=LogisticRegression(),
n_folds=2,
random_state=0,
simultaneous=False,
)
print(result.summary())
Give pseudo_learner a learner of a real-valued target. A transform can leave \([0, 1]\), so a
level can leave the outcome range too. The fit then warns with ParameterSpaceWarning, and
ratio() refuses a scale that is undefined at that level.
The sequentially doubly robust estimator
states the contract and what differs from the default method.
Report the sequential baselines#
SequentialGComputation and SequentialIPW estimate the regimen levels and contrasts that
LTMLE estimates. Report them beside an LTMLE or SDR fit as reference estimators.
estimator |
the estimate reads |
use it to |
|---|---|---|
|
the outcome and pseudo-outcome regressions |
check the targeted estimate against the plug-in |
|
the treatment and censoring mechanisms |
check the targeted estimate against the weighted mean |
from cleverly.datasets import make_longitudinal
from cleverly.learners import GLM
from cleverly.longitudinal import SequentialGComputation, SequentialIPW
frame, _ = make_longitudinal(n=300, seed=0)
columns = dict(
outcome="Y",
treatment=["A1", "A2"],
baseline=["W1", "W2"],
time_varying=[[], ["L2"]],
censoring=["C1", "C2"],
)
regimens = {"always": 1, "never": 0}
gcomp = SequentialGComputation(regimens, reference="never", outcome_learner=GLM()).fit(
frame, **columns
)
ipw = SequentialIPW(regimens, reference="never", treatment_learner=GLM()).fit(frame, **columns)
print(gcomp.summary())
print(ipw.summary())
gcomp() and ipw() are the one-call forms. Both take the same columns and regimens=.
The interval of each baseline holds only when every GLM that the estimate reads is correctly
specified. Under misspecification the interval covers the limit of the misspecified fit, not the
target. A targeted estimator has no such condition on one of its two models.
When a regression that the estimate reads uses a data-adaptive learner, the result reports the
status adaptive_baseline_plugin. It gives the point estimate only. The standard error is a
diagnostic under plugin_std_error. Fit LTMLE or SDR when you need an interval with
data-adaptive learners.
The answers can differ from LTMLE for three reasons.
reason |
effect |
|---|---|
g-computation applies no fluctuation |
its estimate reads the outcome regressions alone |
IPW uses no outcome regression |
its estimate reads the mechanisms alone, and a large weight moves it |
both baselines fit |
|
With the same outcome learner, the g-computation estimate equals the untargeted first step of
LTMLE. Sequential g-computation and IPW gives
the estimators, the stacked influence curve, the settings, and the refusals. The
sequential offers tutorial fits g-computation and LTMLE
with the same misspecified outcome model on one campaign.
Estimate structural blips#
StructuralNestedMean fits a linear structural nested mean model. Use it to report the effect of
each treatment node on the end-of-study outcome, and to let that effect vary with declared
features. It is a separate method. LTMLE and SDR fit no model for the effect.
you want |
fit |
|---|---|
the effect of node \(t\) with later treatment at the reference, as a coefficient matrix |
|
the mean outcome under a static regimen, with no model for the effect |
|
the value of a rule, a policy or an incremental intervention |
|
Declare the features of each node in blip_features. A key is a node counted from one, or its
treatment column. A value lists columns. A tuple of columns declares a product. The intercept is
always a feature. A feature may read the baseline columns, the covariates through its node, and
the earlier treatments.
from cleverly.datasets import make_longitudinal
from cleverly.learners import GLM
from cleverly.longitudinal import StructuralNestedMean
frame, _ = make_longitudinal(n=300, seed=0, censoring=False)
result = StructuralNestedMean(
{"always": 1, "never": 0},
reference="never",
blip_features={1: ["W1"], 2: ["W1", "A1"]},
outcome_learner=GLM(family="binomial"),
treatment_learner=GLM(family="binomial"),
n_folds=2,
random_state=0,
).fit(
frame,
outcome="Y",
treatment=["A1", "A2"],
baseline=["W1", "W2"],
time_varying=[[], ["L2"]],
)
print(result.summary())
print(result.blip(2))
snmm() is the one-call form. result.blip(node) returns the coefficients of a node as a matrix
with one row per feature and one column per coding component. It is the transpose of the matrix
\(\Psi_t\) of the
technical page. A categorical node has one column for each level
after the first sorted level. A continuous node has one column, the dose. The regimen levels and
contrasts follow from the blips and the mean outcome at the reference level ey_reference.
Pass regression_learner= when the outcome is binary and a node holds a dose or a feature is
continuous. A classifier takes no target outside \(\{0, 1\}\), and regression_learner fits those
targets. It defaults to outcome_learner.
Read the result with this condition. The intervals are valid under the linear structural nested
mean model with the declared blip features. Every row of the summary and of to_frame() carries
that label. A blip that depends on an undeclared earlier treatment or covariate makes the
coefficients answer a different question.
The evidence label is “Python property evidence, with comparator runs pending”. The
structural nested mean study
covers cross-fitted fits with GLM nuisances on binary, three-level and dose nodes. It also
covers known weights and 40 or more clusters. A fit inside those rows reads supported. Any other fit, such as an in-sample fit, reads experimental. The summary
prints a support column and the reason.
Check the model against LTMLE.
check = result.check_model()
print(check.summary())
check_model() fits LTMLE on the same rows and folds. A cross-fitted check refits both
estimators on n_splits seeded splits (default 5) and averages the difference. It reports a joint
Wald statistic on the difference of the contrasts. A rejection says that the model fails for these contrasts. A
non-rejection does not confirm the model.
the fit |
what it reports |
|---|---|
|
the influence-curve intervals |
|
the influence-curve intervals, valid when every |
|
the point estimate, with a diagnostic standard error |
The estimator refuses a rule, a policy, censoring, survival and repeated outcomes. It refuses a
static level, a contrast and check_model when a feature reads a time-varying covariate.
Structural dynamic effects lists
each refusal and its alternative.
Simulate the parametric g-formula#
ParametricGFormula estimates the same regimen levels and contrasts as LTMLE. It fits a model for
each time-varying covariate and for the outcome. It then simulates each unit forward under each
regimen and averages the simulated outcomes. SequentialGComputation reaches the same target by
regression backward through the nodes.
reason |
choose |
|---|---|
you need the forward-simulation form that |
|
a plan fixes a continuous dose at a node, or sets it by a known rule |
|
you doubt the law of a time-varying covariate |
|
you want a fit that a misspecified model of one kind cannot bias |
|
The parametric g-formula needs a correct model for every covariate and for the outcome. It reads more models than the backward form, and each one must be right.
from cleverly.datasets import make_longitudinal
from cleverly.longitudinal import ParametricGFormula
frame, _ = make_longitudinal(n=300, seed=0)
result = ParametricGFormula(
{"always": 1, "never": 0},
reference="never",
covariate_types={"L2": "normal"},
n_draws=20,
random_state=0,
).fit(
frame,
outcome="Y",
treatment=["A1", "A2"],
baseline=["W1", "W2"],
time_varying=[[], ["L2"]],
censoring=["C1", "C2"],
)
print(result.summary())
print(result.simulation.monte_carlo_share)
gformula() is the one-call form. The fit refuses an undeclared time-varying column that is not
coded 0/1, before any learner. Declare it in covariate_types= as "binary", "categorical" or
"normal".
The interval includes the simulation error. result.simulation.monte_carlo_share gives the share
of each variance that the simulation contributes, and the fit warns above 10%. Raise n_draws= to
lower it. A fit with the same random_state= repeats exactly.
A data-adaptive learner for any model that the estimate reads gives the status
adaptive_baseline_plugin: the point estimate only. A fixed number at a continuous treatment node
is a fixed dose, and a DynamicRegimen rule with rule_kind="known" there sets the dose from the
simulated history. The result then states its model reliance and reports a support diagnostic. No
published analytic variance exists for the simulation form, so the interval is a natural
extension. Parametric g-formula gives the
estimator, the argument, the settings, the refusals and the validation.
Estimate an incremental intervention over time#
Pass IncrementalEffect to identify() and method="ipsi" to estimate(). Each
Incremental(delta) multiplies every unit’s odds of treatment by delta at every node.
Incremental(1.0) is the natural course, and the example uses it as the reference. A node with
three or more levels takes one factor for each level, such as
Incremental({"high": 2.0, "low": 0.5}). A level that the mapping does not name keeps the
factor 1.
from sklearn.linear_model import LinearRegression, LogisticRegression
from cleverly import CausalStudy, IncrementalEffect, LongitudinalTreatment
from cleverly.datasets import make_longitudinal
from cleverly.interventions import Incremental
frame, _ = make_longitudinal(n=300, seed=0)
study = CausalStudy(
frame,
design=LongitudinalTreatment(
outcome="Y",
treatment=("A1", "A2"),
baseline=("W1", "W2"),
time_varying=((), ("L2",)),
censoring=("C1", "C2"),
),
)
effect = study.identify(IncrementalEffect([Incremental(1.0), Incremental(2.0)]))
result = effect.estimate(
method="ipsi",
outcome_learner=LinearRegression(),
treatment_learner=LogisticRegression(),
censoring_learner=LogisticRegression(),
n_folds=2,
random_state=0,
simultaneous=False,
)
assert result.config.estimator == "ipsi"
assert list(result.estimates) == ["ate_regimen[odds x2 vs natural course]"]
print(result.summary())
The ipsi method fits the one-step estimator of Kennedy (2019), with the dropout weights of Kim,
Kennedy and Naimi (2021). It needs no treatment positivity at a tilt node. It is not doubly
robust, so give treatment_learner a learner that can fit the mechanism at every node.
available_methods() lists tmle and sdr as unavailable for an estimated tilt, with the
reason.
task |
how |
|---|---|
fit the engine directly |
|
use one tilt per node |
|
tilt one node and fix another |
|
tilt a node with three levels |
|
declare a known treatment mechanism |
pass |
A level can leave the outcome range, because the estimate is an average of influence-function
values. The fit then warns with ParameterSpaceWarning. A risk curve can also be non-monotone
across horizons. A mean of a non-absorbing outcome \(Y_t\) is the end-of-study estimand on the panel
truncated after node \(t\). Earlier \(Y_s\) are covariates there (Kim et al. 2021, Theorem 3.1).
A fit with a node of three or more levels carries the label of the registered study
categorical incremental policies
when it matches that study’s rows. Any other such fit is labelled experimental. IPSI is Kennedy’s odds tilt, and lmtp::ipsi is the risk-ratio tilt RiskRatioTilt.
Incremental interventions over time states
the contract.
Longitudinal MSM projections#
MSMProjection can project regimen-specific longitudinal means onto a declared working model. The
regimen grid must identify the coefficient vector; a rank-deficient working design is refused.
A longitudinal MSM projection cross-fits by default, as every longitudinal fit does. Each fold runs an untargeted recursion for every regimen, and one stacked fluctuation per node pools the update. MSM projections states the construction and its evidence.
Diagnostics#
Longitudinal results provide cumulative support, targeting-score, and nuisance-model reports by
node. Call result.diagnostics.support() for leverage and truncation by stage.
The nuisance report separates four roles. Each fitted treatment and censoring model appears once
per node. Outcome and pseudo-outcome models appear per fitted regimen, cause, horizon, and node.
Each row labels its loss as out_of_fold or in_sample.
The report uses weighted negative log likelihood for treatment and censoring. It uses weighted
Brier loss for a binary final target and weighted mean squared error otherwise. Earlier
pseudo-outcomes use weighted mean squared error. A complete-data fit records why no censoring
learner appears in nuisance.omissions.
Use an explicit grid to describe point-estimate movement across cumulative mechanism bounds:
curve = result.diagnostics.truncation_curve(bounds=[0.01, (0.05, 0.95)])
A scalar b means (b, 1), not the point-treatment shorthand (b, 1 - b). The pair bounds each
cumulative treatment and censoring product after multiplication. The diagnostic keeps raw
mechanism predictions fixed and reruns the complete backward recursion at every pair. It refits
every bound-dependent outcome or pseudo-outcome regression and targeting update.
The curve is descriptive. It reports no standard error, confidence interval, preferred bound, or
pass threshold. It refuses mechanism=True, which is a point-treatment option. It also refuses a
learner it cannot replay.
Replay-only unavailability
states the random_state rule each outcome and pseudo-outcome learner must satisfy.
Truncation stability holds the returned frame’s columns and its two score-cell counting rules.
A combined report requires both the grid and permission for expensive refits:
report = result.diagnostics.run_all(
include_refits=True,
arguments={"truncation_curve": {"bounds": [0.01, 0.05]}},
)
Point-only sensitivity formulas remain unavailable unless a longitudinal derivation exists.
Tan (2025) derives population sensitivity bounds for binary, static longitudinal strategies, but no sample estimator for them. The roadmap records that boundary.
Persistence#
result.save() and cleverly.load() carry a longitudinal result through a round trip. The
artifact keeps the folds, the fitted mechanisms, the sequential steps, the targeting state, the
causal metadata, and the truncation replay recipe. It keeps the plan that the fit used. A result
fitted on an adapted panel keeps its PanelLayout and unit_ids, so to_long_panel() works on
the loaded result.
save() writes the result with cloudpickle. A rule, a policy density, a map, or a learner that
you define inline is saved by value. That covers a lambda, a closure, and a function defined in a
notebook or in __main__. A loaded result evaluates each of them again. An importable function
or class is saved by reference, so its module must import on the machine that loads the artifact.
A by-value object loads only under the Python minor version that saved it. Load such an artifact
with the same minor version, for example 3.13 for an artifact that 3.13 wrote.
tests/unit/test_long_panel.py::test_policies_defined_in_main_survive_a_save_in_a_fresh_interpreter
saves a lambda rule, a closure density, and an each_node rule from __main__. It loads them in
a new interpreter and compares the assignments and the density.
Persistence and replayability states the
shared contract for every result, including the assessment cache and the capability rows that a
restored artifact refuses. tests/unit/test_serialization.py’s
test_longitudinal_result_retains_the_complete_fitted_graph_and_assessment checks the round trip
against one cross-fitted, weighted, censored fit.