Time-to-event TMLE: readmission when death competes#

The situation#

A health plan runs a navigation program for adults discharged home from participating hospitals. A navigator can contact a patient at discharge and again on day 31. The plan wants to know if navigation lowers unplanned readmission in the first 60 days after discharge.

Follow-up has two 30-day periods. Three events can end a patient’s follow-up.

event

what it does to the record

role in the analysis

first unplanned readmission

ends the time at risk for the outcome

the event of interest

death before a recorded readmission

makes a later readmission impossible

a competing event

disenrollment from the plan

stops the claims, so the plan cannot see a later readmission or death

censoring

Death is not missing data. A patient who dies cannot be readmitted, and more follow-up does not change that. Disenrollment is different. A disenrolled patient can still be readmitted, but the plan does not see it.

The plan records the variables in this order:

\[ W \rightarrow A_1 \rightarrow C_1 \rightarrow (R_1, D_1) \rightarrow L_2 \rightarrow A_2 \rightarrow C_2 \rightarrow (R_2, D_2) \]

symbol

column

meaning

recorded for

\(W\)

age, baseline_readiness

baseline covariates, measured before discharge

every patient

\(A_1\)

navigation_p1

the navigation offer at discharge

every patient

\(C_1\)

enrolled_p1

1 if the patient stays enrolled through period 1

every patient

\(R_1\), \(D_1\)

readmission_p1, death_p1

readmission, or death before readmission, in period 1

patients enrolled through period 1

\(L_2\)

identified_needs

the care needs recorded by day 30

patients at risk on day 31

\(A_2\)

navigation_p2

the navigation offer on day 31

patients at risk on day 31

\(C_2\), \(R_2\), \(D_2\)

enrolled_p2, readmission_p2, death_p2

the same events in period 2

patients at risk on day 31

Each column is an observed variable. The question is about counterfactual quantities, which are never columns. \(R_t^{a}\) is 1 if the patient would be readmitted by the end of period \(t\) under navigation plan \(a\), had the patient stayed enrolled. \(D_t^{a}\) is the same for death before readmission.

The causal graph#

The first line gives the time order. Each later line lists the direct effects of one node in the synthetic law of this page. The law has these arrows and no others.

time order:  W   A1   C1   (R1, D1)   L2   A2   C2   (R2, D2)

W         -->  A1, C1, R1, D1, L2, A2, R2, D2
A1        -->  C1, R1, D1, L2, A2, R2, D2
L2        -->  A2, C2, R2, D2
A2        -->  R2, D2
(R1, D1)  -->  every period-2 node: an event ends follow-up
C1, C2    -->  the record: disenrollment hides the later events

arrow or absent arrow

what it asserts

consequence for the analysis

\(W \rightarrow A_1, A_2\)

age and readiness change who receives an offer

\(W\) confounds both offers

\(W \rightarrow R_t, D_t, C_1, L_2\)

the same covariates change the event hazards. age alone also changes period-1 disenrollment and the needs

the analysis adjusts for \(W\) at every node

\(A_1 \rightarrow L_2\)

navigation at discharge raises the recorded needs

\(L_2\) is a mediator of the first offer

\(L_2 \rightarrow A_2, R_2, D_2\)

needs drive the day-31 offer and the period-2 hazards

\(L_2\) confounds the second offer

\(A_1 \rightarrow C_1\), \(L_2 \rightarrow C_2\)

navigation and needs change disenrollment

censoring depends on the history, so the analysis must model it

\(A_1 \rightarrow R_t, D_t\), \(A_2 \rightarrow R_2, D_2\)

navigation changes both hazards and which cause occurs

these paths make up the effect to estimate

\((R_1, D_1) \rightarrow\) period 2

an event in period 1 ends the time at risk

readmission and death are absorbing, and death competes with readmission

no unmeasured common cause of an offer, disenrollment, or an event

the recorded history holds every common cause

this absence is sequential exchangeability and independent censoring. No data check can confirm it

no arrow between patients

one patient’s offer does not change another patient’s events

this absence is no interference

\(L_2\) is a mediator for the first offer and a confounder for the second. One regression on both offers cannot handle it. Adjustment for \(L_2\) blocks part of the effect of the first offer. Omission of \(L_2\) leaves the second offer confounded. A sequential method adjusts for \(L_2\) at the second decision only.

The hazards depend on \(W\) and \(L_2\), so the effect can differ between patients. These variables are effect modifiers. The estimands below average over them.

The question and the estimand#

Target population. Adults discharged home from a participating hospital and enrolled in the plan at discharge.

Interventions. The analysis compares two plans.

plan

rule

always

offer navigation at discharge, and again on day 31 if the patient has had neither event

never

offer no navigation at either decision

Under both plans, the estimand describes the patients had all of them stayed enrolled. This is a hypothetical strategy for censoring. It suits disenrollment, because disenrollment hides events and does not prevent them.

Estimands. For plan \(a\) and period \(t\), where \(t = 1\) ends on day 30 and \(t = 2\) on day 60:

estimand

definition

what it answers

readmission incidence \(F_R^{a}(t)\)

\(P(R_t^{a} = 1)\): readmission before death, by period \(t\)

how many patients the plan would see readmitted

death incidence \(F_D^{a}(t)\)

\(P(D_t^{a} = 1)\): death before readmission, by period \(t\)

how many patients die first

event-free survival \(S^{a}(t)\)

\(1 - F_R^{a}(t) - F_D^{a}(t)\)

how many patients reach the end of period \(t\) alive and not readmitted

The effects are the always-minus-never differences. The question is about \(F_R\).

The analysis also reports \(F_D\) and \(S\), for one reason. If navigation prevents deaths, more patients remain at risk of readmission. \(F_R\) can then rise even when navigation lowers each patient’s readmission hazard. Young et al. (2020) call the readmission difference a total effect. It includes this path through death.

Why a simpler estimand does not answer the question.

simpler estimand

why it fails here

where the page shows it

the readmission share among patients whose offers followed a plan

offers go to patients at higher risk, and disenrolled patients drop out

Step 3

readmission risk with death treated as censoring, such as one minus Kaplan-Meier

it describes a world in which nobody dies. That world needs more assumptions, and the program cannot create it

Step 8

the effect of the discharge offer alone

it describes a different intervention, and \(L_2\) makes one regression wrong for two offers

the causal graph

Identification assumptions.

assumption

what it requires here

can the data check it?

consistency

each recorded offer is the declared navigation version, and each recorded event is the event under that plan

no. It rests on the protocol’s version records

no interference

one patient’s offers do not change another patient’s events

no. It rests on the navigator capacity design

sequential exchangeability

at each decision, among patients still event-free and enrolled, the offer is independent of the counterfactual events, given \(W\), \(L_2\), and earlier offers

no. It rests on judgment about the recorded history

independent censoring

in each period, staying enrolled is independent of the counterfactual events, given the same history

no. It rests on the same judgment, for disenrollment

positivity

every kind of patient has some chance of following each plan and staying enrolled

partly, through the support report of Step 9

Under these assumptions, the sequential g-formula identifies each estimand from the observed data.

Why this method#

Each feature of the method answers one feature of the problem.

feature of the problem

method feature

what it costs

two offers, and \(L_2\) is caused by the first and confounds the second

longitudinal TMLE: one regression at each node, from the last period back to the first

one regression per node, per horizon, per cause

readmission and death are absorbing

each node’s regression uses only the patients at risk entering that node

none

death competes with readmission

a cause-specific recursion. A patient who died carries a zero forward for readmission

one backward pass per cause. The causes share the treatment and censoring models

disenrollment depends on the history

the targeting step weights by the probability of following the plan and staying enrolled

it needs positivity for staying enrolled

4,000 patients and compact parametric learners

in-sample fits without cross-fitting

inference rests on the Donsker condition, which these learners satisfy

For one event type, outcome= takes one column per period. For competing events, it takes a mapping from each cause to its columns.

term

plain meaning

node

one time at which a column is recorded, such as the day-31 offer

risk set

the patients who have had neither event and are still enrolled when a node starts

nuisance

a model the estimate needs but the question does not ask about, such as the offer model

targeting

a small update to each outcome regression, weighted by the treatment and censoring models, that removes first-order bias

influence curve

how much each patient moves the estimate. Its variance gives the standard error

Targeting and bounds and the inference entry define the last two terms. The technical entry gives the event-process recursion.

Step 1: set up#

The setup imports the learners and prints the installed cleverly version. Every fit passes its learners and random seed explicitly, so a rerun reproduces the stored outputs.

import pandas as pd
from sklearn.linear_model import LinearRegression, LogisticRegression

import cleverly

pd.set_option("display.width", 160)
pd.set_option("display.max_columns", 20)
print("cleverly", cleverly.__version__)
cleverly 0.1.3

What this output tells you. The stored outputs came from the version named above. A different version can print different numbers.

Step 2: the data#

The data come from make_longitudinal_competing, a synthetic law with a known answer. The code renames its columns to the program’s names. It counts how patients leave the first risk set, and it prints the true value of each estimand.

from cleverly.datasets import make_longitudinal_competing

frame, truth = make_longitudinal_competing(n=4_000, seed=53, censoring=True)
frame = frame.rename(
    columns={
        "W1": "age",
        "W2": "baseline_readiness",
        "A1": "navigation_p1",
        "C1": "enrolled_p1",
        "R1": "readmission_p1",
        "D1": "death_p1",
        "L2": "identified_needs",
        "A2": "navigation_p2",
        "C2": "enrolled_p2",
        "R2": "readmission_p2",
        "D2": "death_p2",
    }
)
print("rows and columns:", frame.shape)
print(list(frame.columns))
print(frame.head().round(3).to_string())
print()
print("disenrolled in period 1:   ", int(frame["enrolled_p1"].eq(0).sum()))
print("readmitted in period 1:    ", int(frame["readmission_p1"].eq(1).sum()))
print("died in period 1:          ", int(frame["death_p1"].eq(1).sum()))
print("at risk entering period 2: ", int(frame["navigation_p2"].notna().sum()))
print("disenrolled in period 2:   ", int(frame["enrolled_p2"].eq(0).sum()))
print()

# The generator names readmission "relapse". The page uses the program's word.
population = {key.replace("relapse", "readmission"): value for key, value in truth.items()}
rows = []
for horizon in (1, 2):
    row = {"time": horizon}
    for plan in ("always", "never"):
        readmission = population[f"cif_regimen[{plan}, readmission @ t={horizon}]"]
        death = population[f"cif_regimen[{plan}, death @ t={horizon}]"]
        row |= {
            f"readmission_{plan}": readmission,
            f"death_{plan}": death,
            f"event_free_{plan}": 1.0 - readmission - death,
        }
    rows.append(row)
print("known values of the synthetic law, by plan:")
print(pd.DataFrame(rows).round(4).to_string(index=False))
rows and columns: (4000, 11)
['age', 'baseline_readiness', 'navigation_p1', 'enrolled_p1', 'readmission_p1', 'death_p1', 'identified_needs', 'navigation_p2', 'enrolled_p2', 'readmission_p2', 'death_p2']
     age  baseline_readiness  navigation_p1  enrolled_p1  readmission_p1  death_p1  identified_needs  navigation_p2  enrolled_p2  readmission_p2  death_p2
0  0.105               0.891            0.0          1.0             0.0       0.0            -0.347            0.0          1.0             1.0       0.0
1 -0.807              -0.103            1.0          0.0             NaN       NaN               NaN            NaN          NaN             NaN       NaN
2  0.892               1.262            0.0          1.0             0.0       0.0             0.102            1.0          1.0             0.0       0.0
3 -0.929               0.727            0.0          1.0             0.0       0.0            -1.423            1.0          1.0             0.0       0.0
4 -1.136              -0.282            0.0          1.0             1.0       0.0               NaN            NaN          NaN             1.0       0.0

disenrolled in period 1:    451
readmitted in period 1:     457
died in period 1:           253
at risk entering period 2:  2839
disenrolled in period 2:    233

known values of the synthetic law, by plan:
 time  readmission_always  death_always  event_free_always  readmission_never  death_never  event_free_never
    1              0.1187        0.0310             0.8503             0.1452       0.1127            0.7421
    2              0.2467        0.0705             0.6828             0.2443       0.2109            0.5448

What this output tells you. Each row is one patient. Row 1 disenrolled in period 1, so every later node is NaN. Row 4 was readmitted in period 1, so readmission_p2 carries the 1 forward and the day-31 offer does not exist.

In period 1, 451 patients disenrolled, 457 were readmitted, and 253 died. That leaves 2839 patients at risk entering period 2. Of these, 233 disenrolled in period 2.

plan

readmission by day 60

death by day 60

event-free at day 60

always

0.2467

0.0705

0.6828

never

0.2443

0.2109

0.5448

In this law, navigation lowers the hazard of any event and raises the share of events that are readmissions. Death by day 60 falls from 0.2109 to 0.0705. Readmission by day 30 falls from 0.1452 to 0.1187, but readmission by day 60 rises from 0.2443 to 0.2467. A real program has no known values. The comparisons with them below are a teaching device.

Step 3: association first#

This comparison shows why the analysis needs adjustment. Its target is the always-minus-never difference of each incidence. The naive estimate is the event share among patients whose recorded offers followed each plan and who stayed enrolled. It assumes that the offers and disenrollment are unrelated to risk. The causal graph says that they are related, through \(W\) and \(L_2\).

print("baseline means by the discharge offer:")
print(frame.groupby("navigation_p1")[["age", "baseline_readiness"]].mean().round(3))
print()
at_risk_p2 = frame["navigation_p2"].notna()
needs = frame.loc[at_risk_p2].groupby("navigation_p2")["identified_needs"].mean()
print("mean identified_needs by the day-31 offer, among patients at risk:")
print(needs.round(3))
print()

# A follower's recorded offers match the plan at every decision the patient reached.
first = {"always": frame["navigation_p1"].eq(1), "never": frame["navigation_p1"].eq(0)}
second = {"always": ~frame["navigation_p2"].eq(0), "never": ~frame["navigation_p2"].eq(1)}
followers = {1: first, 2: {plan: first[plan] & second[plan] for plan in first}}
naive_differences = {}
print("share with the event among enrolled followers, always minus never:")
for cause in ("readmission", "death"):
    for horizon in (1, 2):
        outcome = frame[f"{cause}_p{horizon}"]
        shares = {
            plan: float(outcome[rows & outcome.notna()].mean())
            for plan, rows in followers[horizon].items()
        }
        naive = naive_differences[cause, horizon] = shares["always"] - shares["never"]
        target = population[f"ate_regimen[always vs never, {cause} @ t={horizon}]"]
        print(
            f"  {cause:11s} t={horizon}: always {shares['always']:.4f}, "
            f"never {shares['never']:.4f}, difference {naive:.4f}, population {target:.4f}"
        )
baseline means by the discharge offer:
                 age  baseline_readiness
navigation_p1                           
0.0           -0.127               0.179
1.0            0.140              -0.197

mean identified_needs by the day-31 offer, among patients at risk:
navigation_p2
0.0   -0.104
1.0    0.793
Name: identified_needs, dtype: float64

share with the event among enrolled followers, always minus never:
  readmission t=1: always 0.1345, never 0.1236, difference 0.0109, population -0.0265
  readmission t=2: always 0.3317, never 0.2707, difference 0.0610, population 0.0024
  death       t=1: always 0.0314, never 0.1075, difference -0.0761, population -0.0817
  death       t=2: always 0.0920, never 0.2425, difference -0.1506, population -0.1404

What this output tells you. The groups differ before each offer.

comparison

offered

not offered

mean age at the discharge offer

0.140

-0.127

mean baseline_readiness at the discharge offer

-0.197

0.179

mean identified_needs at the day-31 offer

0.793

-0.104

In this law, older age, lower readiness, and more needs raise the event hazard. The offers therefore go to patients at higher risk.

cause and day

naive difference

population difference

readmission, day 30

0.0109

-0.0265

readmission, day 60

0.0610

0.0024

death, day 30

-0.0761

-0.0817

death, day 60

-0.1506

-0.1404

On this draw, the naive day-30 readmission difference has the wrong sign. The naive share also drops every disenrolled patient. Disenrollment depends on navigation and needs, so the patients who remain are not a random subset.

Step 4: write the protocol#

A StudyProtocol records the scientific design before any model runs. Point-treatment TMLE introduces its fields. The code starts from longitudinal_navigation_protocol(), the two-decision program of the shared study design. It replaces only the fields that this question changes.

from dataclasses import replace

from cleverly.datasets import longitudinal_navigation_protocol

program = longitudinal_navigation_protocol()
protocol = replace(
    program,
    eligibility=(*program.eligibility, "Enrolled in the health plan at discharge"),
    treatment_strategies=(
        "Offer navigation at the start of each 30-day period until readmission or death",
        program.treatment_strategies[1],
    ),
    treatment_versions=(
        "The declared navigation contacts at discharge and on day 31",
        program.treatment_versions[1],
    ),
    outcome=(
        "First unplanned readmission, and death before a recorded readmission, "
        "recorded at the end of each 30-day period"
    ),
    horizon="30 and 60 days after discharge",
    intercurrent_event_handling=(
        "Treat death before readmission as a competing event that stays in the history",
        "Treat disenrollment before readmission or death as censoring, so the target is the "
        "incidence had every patient stayed enrolled (hypothetical strategy)",
        "Make no navigation decision after either event",
        program.intercurrent_event_handling[2],
    ),
    assumption_rationale=(
        "Recorded history covers the measured common causes of navigation, disenrollment, "
        "readmission, and death in each period",
        *program.assumption_rationale[1:],
    ),
)
print("\n".join(protocol.summary_lines()))
causal study protocol: schema 1; 9b9cb1dda3683cf0
target population: Adults discharged home from a participating hospital during the enrollment period
eligibility: ['Age 18 years or older', 'Discharged alive', 'Discharged home from a participating hospital', 'Enrolled in the health plan at discharge']
time zero: Hospital discharge, after baseline measurement and before first assignment
treatment strategies: ['Offer navigation at the start of each 30-day period until readmission or death', 'Offer no navigation at either decision']
treatment versions: ['The declared navigation contacts at discharge and on day 31', 'Usual discharge support without navigation contacts']
outcome: First unplanned readmission, and death before a recorded readmission, recorded at the end of each 30-day period
horizon: 30 and 60 days after discharge
intercurrent-event handling: ['Treat death before readmission as a competing event that stays in the history', 'Treat disenrollment before readmission or death as censoring, so the target is the incidence had every patient stayed enrolled (hypothetical strategy)', 'Make no navigation decision after either event', 'Analyze each navigation offer regardless of completed contacts']
interference unit: Individual patient
assumption rationale: ['Recorded history covers the measured common causes of navigation, disenrollment, readmission, and death in each period', 'Version records support consistency at both navigation decisions', 'Reserved navigator capacity supports no interference between patients']

What this output tells you. The fingerprint is 9b9cb1dda3683cf0. The code changes seven fields: the eligibility, the treatment strategies and versions, the outcome, the horizon, the intercurrent-event handling, and the assumption rationale.

The intercurrent-event handling records the decisions of the estimand section. Death before readmission is a competing event. Disenrollment is censoring under a hypothetical strategy. No navigation decision follows either event. The protocol stores text only, so the design of Step 5 must match these entries.

Step 5: design and identification#

The design places every column at its node. The outcome= mapping declares two competing causes, with one indicator column per period for each. The censoring= role declares enrollment in each period. time_varying= places identified_needs before the day-31 offer. RegimeMean asks for each incidence under each plan at both horizons. The estimands guide defines the estimand objects.

from cleverly import CausalStudy, LongitudinalTreatment, RegimeMean

# Step 7 reuses the offers and the history with another censoring definition.
navigation_nodes = {
    "treatment": ("navigation_p1", "navigation_p2"),
    "baseline": ("age", "baseline_readiness"),
    "time_varying": ((), ("identified_needs",)),
}
plans = {"always": 1, "never": 0}
levels_by_plan = RegimeMean(plans, reference="never", horizons=(1, 2))
study = CausalStudy(
    frame,
    design=LongitudinalTreatment(
        outcome={
            "readmission": ("readmission_p1", "readmission_p2"),
            "death": ("death_p1", "death_p2"),
        },
        censoring=("enrolled_p1", "enrolled_p2"),
        **navigation_nodes,
    ),
    protocol=protocol,
)
effect = study.identify(levels_by_plan)
print(effect.summary(protocol="fingerprint"))
mean outcome under each declared regime
identified by explicit-adjustment: sequential g-formula under the declared treatment regimen
history at navigation_p1: ['age', 'baseline_readiness']
history at navigation_p2: ['age', 'baseline_readiness', 'identified_needs', 'navigation_p1']
required nuisances: ['sequential_outcome_regressions', 'treatment_and_censoring_mechanisms']
assumptions:
  - consistency: each observed history equals the potential history under its realized regimen
  - no interference: one unit's potential history does not depend on other units' regimens
  - sequential exchangeability given the recorded history at every node
  - sequential positivity for treatment and remaining under observation
causal study protocol: schema 1; 9b9cb1dda3683cf0

What this output tells you. The identification is the sequential g-formula. The history at the day-31 offer is age, baseline_readiness, identified_needs, and the first offer. The first offer has no identified_needs, because the needs are recorded later. The required nuisances are the sequential outcome regressions and the treatment-and-censoring mechanisms. The printed assumptions are those of the estimand section. “Remaining under observation” is the independent-censoring assumption for enrollment.

Step 6: estimate the incidences#

The fit uses compact parametric learners on all rows, without cross-fitting.

setting

value

why

outcome_learner

logistic regression

the last node regresses a binary event

pseudo_learner

linear regression

the earlier node regresses a predicted probability, not an event

treatment_learner, censoring_learner

logistic regression on the history

in this law, these have the form of the true offer and enrollment models

CrossFitting(enabled=False)

in-sample

these learners satisfy the Donsker condition, so in-sample fits keep influence-curve inference valid under the remaining conditions. With misspecified outcome regressions, the standard error can be conservative

Runtime(random_state=41, n_jobs=1)

fixed seed, one process

the output is reproducible

The outcome regressions do not match the law. Its hazards contain a tanh term, and each cause is a hazard times a cause share. TMLE is doubly robust: the estimate stays consistent when the offer and enrollment models are correct. Here that condition holds by construction. A real analysis cannot verify it.

from cleverly import CrossFitting, ModelSpec, Runtime, TMLEMethod

sequential = TMLEMethod(
    models=ModelSpec(
        outcome_learner=LogisticRegression(max_iter=1000, random_state=41),
        pseudo_learner=LinearRegression(),
        treatment_learner=LogisticRegression(max_iter=1000, random_state=41),
        censoring_learner=LogisticRegression(max_iter=1000, random_state=41),
    ),
    cross_fitting=CrossFitting(enabled=False),
    runtime=Runtime(random_state=41, n_jobs=1),
)
result = effect.estimate(method=sequential)
print(result.summary(protocol="fingerprint"))
Longitudinal TMLE (2 time points, n = 4000)

parameter                               estimate  std. error  95% CI            p-value
--------------------------------------  --------  ----------  ----------------  -------
cif_regimen[always, readmission @ t=1]  0.1192    0.0076      [0.1042, 0.1341]  <1e-4  
cif_regimen[always, readmission @ t=2]  0.2508    0.0116      [0.2281, 0.2735]  <1e-4  
cif_regimen[always, death @ t=1]        0.0303    0.0042      [0.0220, 0.0385]  <1e-4  
cif_regimen[always, death @ t=2]        0.0722    0.0070      [0.0585, 0.0859]  <1e-4  
cif_regimen[never, readmission @ t=1]   0.1364    0.0084      [0.1200, 0.1529]  <1e-4  
cif_regimen[never, readmission @ t=2]   0.2285    0.0141      [0.2008, 0.2561]  <1e-4  
cif_regimen[never, death @ t=1]         0.1113    0.0076      [0.0965, 0.1262]  <1e-4  
cif_regimen[never, death @ t=2]         0.2098    0.0140      [0.1824, 0.2373]  <1e-4  

  time points: 2
  outcome: competing risks, absorbing causes readmission, death at readmission_p1, readmission_p2, death_p1, death_p2
  horizons reported: t = 1, 2
  causes reported: readmission, death
  outcome family: binomial
  regimens: always=(1/1), never=(0/0)
  cross-fitting: none -- nuisances fitted in sample
  g_bounds: fixed [0.01, 1] on each cumulative treatment-and-censoring probability (package default; R ltmle-compatible heuristic -- inspect truncation diagnostics)
  survival model: iterated_mean
  targeting: iterative
  isotonic projection: none
  confidence level: 95%
  random_state: 41
  causal estimand: mean outcome under each declared regime
  identification: explicit-adjustment; sequential g-formula under the declared treatment regimen
  required nuisances: sequential_outcome_regressions, treatment_and_censoring_mechanisms
  identification assumptions: consistency: each observed history equals the potential history under its realized regimen; no interference: one unit's potential history does not depend on other units' regimens; sequential exchangeability given the recorded history at every node; sequential positivity for treatment and remaining under observation
  causal study protocol: schema 1; 9b9cb1dda3683cf0
  data 0243e934597a0cc0 | folds 4dc224c385945072 | cleverly 0.1.3 | 2026-10-16T02:37:47+00:00
  always: 938 of 4000 units followed it throughout; max weight 14.5, effective n 786; max truncated share 0.0% at t=1
  never: 671 of 4000 units followed it throughout; max weight 28.7, effective n 529; max truncated share 0.0% at t=1

  simultaneous 95% bands (multiplier bootstrap, critical value 2.679 vs 1.960 pointwise):
    cif_regimen[always, readmission @ t=1]  [0.0987, 0.1396]
    cif_regimen[always, readmission @ t=2]  [0.2197, 0.2818]
    cif_regimen[always, death @ t=1]  [0.0190, 0.0415]
    cif_regimen[always, death @ t=2]  [0.0535, 0.0910]
    cif_regimen[never, readmission @ t=1]  [0.1139, 0.1589]
    cif_regimen[never, readmission @ t=2]  [0.1907, 0.2663]
    cif_regimen[never, death @ t=1]  [0.0911, 0.1316]
    cif_regimen[never, death @ t=2]  [0.1723, 0.2473]

What this output tells you. The table has one row per plan, per cause, per horizon. Each row is a cumulative incidence with its pointwise 95% interval. These lines of the summary matter for this fit.

summary line

what it says

causes reported: readmission, death

one backward pass for each cause

cross-fitting: none

the in-sample construction

g_bounds: fixed [0.01, 1]

the lower bound on each cumulative probability of following the plan and staying enrolled. Step 9 shows that no patient reached it

always: 938 of 4000 units followed it throughout

938 patients followed always, and 671 followed never. The largest weights are 14.5 and 28.7

The summary also prints simultaneous bands. This page reads only the pointwise intervals. Method evidence links the registered studies, which also measure the bands. The question asks for differences, which Step 7 computes.

Step 7: compare the plans#

Each difference uses result.contrast(...). It applies the delta method to the joint influence curve of the estimates it combines. All the incidences come from one sample, and the two incidences of one plan come from the same patients. The joint curve includes the correlation between their errors. Event-free survival is one minus the sum of both incidences. One minus a single incidence is not event-free survival.

def incidence(plan, cause, horizon):
    return f"cif_regimen[{plan}, {cause} @ t={horizon}]"


def report(label, estimate, target):
    low, high = estimate.ci
    print(
        f"  {label:22s} {estimate.psi:7.4f}  SE {estimate.std_error:.4f}  "
        f"CI ({low:7.4f}, {high:7.4f})  population {target:7.4f}"
    )


# Each contrast reads the joint influence curve of the estimates it combines.
differences = {}
print("always minus never, cause-specific cumulative incidence:")
for cause in ("readmission", "death"):
    for horizon in (1, 2):
        names = [incidence("always", cause, horizon), incidence("never", cause, horizon)]
        differences[cause, horizon] = result.contrast(
            lambda psi: psi[0] - psi[1], names, name=f"{cause} difference @ t={horizon}"
        )
        target = population[f"ate_regimen[always vs never, {cause} @ t={horizon}]"]
        report(f"{cause} t={horizon}", differences[cause, horizon], target)
print()

# Event-free survival is one minus the sum of both causes, never one minus one cause.
event_free = {}
print("event-free survival, one minus both incidences:")
for horizon in (1, 2):
    for plan in ("always", "never"):
        names = [incidence(plan, "readmission", horizon), incidence(plan, "death", horizon)]
        event_free[plan, horizon] = result.contrast(
            lambda psi: 1.0 - psi[0] - psi[1], names, name=f"event-free {plan} @ t={horizon}"
        )
        target = 1.0 - sum(population[name] for name in names)
        report(f"{plan} t={horizon}", event_free[plan, horizon], target)
    names = [
        incidence(plan, cause, horizon)
        for plan in ("always", "never")
        for cause in ("readmission", "death")
    ]
    differences["event-free", horizon] = result.contrast(
        lambda psi: (psi[2] + psi[3]) - (psi[0] + psi[1]),
        names,
        name=f"event-free difference @ t={horizon}",
    )
    target = sum(population[name] for name in names[2:]) - sum(
        population[name] for name in names[:2]
    )
    report(f"difference t={horizon}", differences["event-free", horizon], target)
always minus never, cause-specific cumulative incidence:
  readmission t=1        -0.0172  SE 0.0113  CI (-0.0393,  0.0048)  population -0.0265
  readmission t=2         0.0223  SE 0.0181  CI (-0.0132,  0.0578)  population  0.0024
  death t=1              -0.0811  SE 0.0087  CI (-0.0980, -0.0641)  population -0.0817
  death t=2              -0.1376  SE 0.0156  CI (-0.1682, -0.1070)  population -0.1404

event-free survival, one minus both incidences:
  always t=1              0.8506  SE 0.0085  CI ( 0.8339,  0.8672)  population  0.8503
  never t=1               0.7522  SE 0.0103  CI ( 0.7320,  0.7725)  population  0.7421
  difference t=1          0.0983  SE 0.0133  CI ( 0.0723,  0.1243)  population  0.1082
  always t=2              0.6770  SE 0.0126  CI ( 0.6523,  0.7017)  population  0.6828
  never t=2               0.5617  SE 0.0167  CI ( 0.5290,  0.5944)  population  0.5448
  difference t=2          0.1153  SE 0.0207  CI ( 0.0748,  0.1558)  population  0.1380

What this output tells you.

always minus never

day 30 (95% CI)

day 60 (95% CI)

population, day 60

readmission incidence

-0.0172 (-0.0393, 0.0048)

0.0223 (-0.0132, 0.0578)

0.0024

death incidence

-0.0811 (-0.0980, -0.0641)

-0.1376 (-0.1682, -0.1070)

-0.1404

event-free survival

0.0983 (0.0723, 0.1243)

0.1153 (0.0748, 0.1558)

0.1380

The readmission intervals contain zero at both horizons. The death and event-free intervals exclude zero. On this draw, each of the ten printed intervals contains its population value. That is one draw, not a coverage result. Method evidence links the repeated-sampling evidence.

Step 8: the failure mode, coding death as censoring#

A common shortcut codes death as censoring, as one minus Kaplan-Meier does. cleverly refuses to eliminate a competing event, with the code different_estimand. A manual recoding of the columns bypasses that refusal. The fit then runs, but it targets a different estimand: readmission risk in a world where death is prevented. The two analyses should differ here for a known reason. Censoring at death treats each patient who died as if the patient stayed at risk of readmission. More patients die under the never plan, so the shortcut raises the never-plan risk the most.

The code builds a new censoring indicator, \(K_t = C_t (1 - D_t)\). A patient stays in the readmission risk set only while enrolled and alive. The code fits the shortcut and writes the protocol it implies. The synthetic law defines no world without death, because it draws an event and then its cause. The code therefore writes the value of the g-formula functional that the recoding targets. It computes that value by quadrature from the generator’s structural equations.

# A patient leaves the readmission risk set at disenrollment or at death, whichever is first.
# The new censoring indicator is K_k = C_k (1 - D_k).
kept_p1 = frame["enrolled_p1"] * (1.0 - frame["death_p1"].fillna(0.0))
through_p1 = kept_p1.eq(1) & frame["readmission_p1"].eq(0)
kept_p2 = (frame["enrolled_p2"] * (1.0 - frame["death_p2"].fillna(0.0))).where(through_p1)
death_as_censoring = frame.assign(
    enrolled_alive_p1=kept_p1,
    enrolled_alive_p2=kept_p2,
    readmission_p1=frame["readmission_p1"].where(kept_p1.eq(1)),
    readmission_p2=frame["readmission_p2"].where(kept_p2.eq(1) | frame["readmission_p1"].eq(1)),
)
eliminated_protocol = replace(
    protocol,
    outcome="First unplanned readmission, recorded at the end of each 30-day period",
    intercurrent_event_handling=(
        "Treat death as censoring, so the target is readmission risk if death were "
        "prevented (hypothetical strategy)",
        protocol.intercurrent_event_handling[1],
        "Make no navigation decision after readmission or death",
        protocol.intercurrent_event_handling[3],
    ),
    assumption_rationale=(
        "Recorded history covers the measured common causes of navigation, disenrollment, "
        "readmission, and death in each period, including every common cause of death and "
        "readmission",
        *protocol.assumption_rationale[1:],
    ),
)
eliminated = (
    CausalStudy(
        death_as_censoring,
        design=LongitudinalTreatment(
            outcome=("readmission_p1", "readmission_p2"),
            censoring=("enrolled_alive_p1", "enrolled_alive_p2"),
            **navigation_nodes,
        ),
        protocol=eliminated_protocol,
    )
    .identify(levels_by_plan)
    .estimate(method=sequential)
)
eliminated_t2 = eliminated.contrast(
    lambda psi: psi[0] - psi[1],
    ["risk_regimen[always @ t=2]", "risk_regimen[never @ t=2]"],
    name="death-as-censoring difference @ t=2",
)

# The law defines no world without death. These constants are values of the g-formula functional
# that the recoding targets, by quadrature over the generator's structural equations. Within a
# period, the functional needs an order for death and readmission, and the two orders differ.
functional = {
    "death first": {"never": 0.304906, "always": 0.264627},
    "readmission first": {"never": 0.260970, "always": 0.251680},
}
print("protocol fingerprint, death as a competing event:", protocol.fingerprint)
print("protocol fingerprint, death as censoring:        ", eliminated_protocol.fingerprint)
print()
print("readmission risk at t=2      never   always  always - never  SE      target")
rows = (
    (
        "death as a competing event",
        result[incidence("never", "readmission", 2)].psi,
        result[incidence("always", "readmission", 2)].psi,
        differences["readmission", 2],
        population["ate_regimen[always vs never, readmission @ t=2]"],
    ),
    (
        "death as censoring",
        eliminated["risk_regimen[never @ t=2]"].psi,
        eliminated["risk_regimen[always @ t=2]"].psi,
        eliminated_t2,
        functional["death first"]["always"] - functional["death first"]["never"],
    ),
)
for label, never, always, contrast, target in rows:
    print(
        f"{label:28s} {never:.4f}  {always:.4f}  {contrast.psi:7.4f}         "
        f"{contrast.std_error:.4f}  {target:7.4f}"
    )
print()
print("death as censoring, 95% CI of the difference:", tuple(round(v, 4) for v in eliminated_t2.ci))
for order, values in functional.items():
    print(f"functional difference, {order:17s}: {values['always'] - values['never']:.4f}")
protocol fingerprint, death as a competing event: 9b9cb1dda3683cf0
protocol fingerprint, death as censoring:         cd5958d44e088c86

readmission risk at t=2      never   always  always - never  SE      target
death as a competing event   0.2285  0.2508   0.0223         0.0181   0.0024
death as censoring           0.2879  0.2698  -0.0180         0.0210  -0.0403

death as censoring, 95% CI of the difference: (-0.0592, 0.0231)
functional difference, death first      : -0.0403
functional difference, readmission first: -0.0093

What this output tells you. The two analyses carry different protocol fingerprints, 9b9cb1dda3683cf0 and cd5958d44e088c86. The death-as-censoring protocol changes three fields: the outcome, the intercurrent-event handling, and the assumption rationale. Both designs run. Only the protocol records that the question changed.

analysis

estimand

day-60 difference (SE)

value of its target

added assumptions

death as a competing event

the total effect on readmission incidence

0.0223 (0.0181)

0.0024

none beyond Step 5

death as censoring

readmission risk if death were prevented

-0.0180 (0.0210)

-0.0403

exchangeability and positivity for death given the history, and a defined intervention that prevents death

On this draw, the shortcut reports a reduction, and the competing-event analysis reports an increase. The never-plan risk rises from 0.2285 to 0.2879. The always-plan risk rises from 0.2508 to 0.2698. The shortcut’s interval, (-0.0592, 0.0231), contains both target values, 0.0024 and -0.0403. This draw cannot tell them apart. The protocol states which estimand the analysis answers.

The shortcut is not doubly robust here. Its censoring indicator mixes disenrollment with death, so the logistic model of \(K_t\) is misspecified, and so are the outcome regressions. Nothing ensures that the shortcut is consistent for its own functional.

Do not read -0.0180 as a direct effect of navigation on readmission. That reading needs a world where death is prevented, and the program has no such intervention.

In discrete time, the functional also depends on the order of death and readmission within a period. With death first, as the recoding assumes, its value is -0.0403. With readmission first, it is -0.0093. The two orders encode two different assumptions about the prevented deaths. A real analysis must defend its order, and no data check can confirm it.

No registered study covers this construction. Young et al. (2020) define both estimands and their identifying conditions.

Step 9: diagnostics, what the fit can check#

The assessment lists the checks that ran on the fit. The support report describes the weights at each node. It counts the patients who follow each plan, gives the largest inverse-probability weight, and says whether any probability reached the bound.

assessment = result.assess()
ledger = assessment.to_frame()
ran = ledger[ledger["status"].isin(["passed", "completed"])]
print(ran[["surface", "check", "status"]].to_string(index=False))
print("needs attention:", tuple(item.name for item in assessment.attention))
print()
support = assessment.report("support").to_frame()
print(support.round(4).to_string(index=False))
ratio = (support["effective_n"] / support["n_followed"]).min()
print(f"minimum effective-sample-size ratio: {ratio:.1%}")
   surface           check    status
validation score_equations    passed
validation         support completed
validation nuisance_models completed
needs attention: ()

regimen       cause  horizon  time  n_followed  share_assigned_1  max_weight  effective_n  share_truncated  epsilon  converged
 always readmission        1     1        1688               1.0      7.0535    1564.5728              0.0  -0.0041       True
 always readmission        2     1        1688               1.0      7.0535    1564.5728              0.0   0.0029       True
 always readmission        2     2         938               1.0     14.5333     786.2736              0.0  -0.0049       True
 always       death        1     1        1688               1.0      7.0535    1564.5728              0.0  -0.0126       True
 always       death        2     1        1688               1.0      7.0535    1564.5728              0.0  -0.0016       True
 always       death        2     2         938               1.0     14.5333     786.2736              0.0  -0.0090       True
  never readmission        1     1        1861               0.0      6.1974    1764.7557              0.0   0.0054       True
  never readmission        2     1        1861               0.0      6.1974    1764.7557              0.0   0.0054       True
  never readmission        2     2         671               0.0     28.7198     529.0079              0.0  -0.0679       True
  never       death        1     1        1861               0.0      6.1974    1764.7557              0.0  -0.0053       True
  never       death        2     1        1861               0.0      6.1974    1764.7557              0.0  -0.0028       True
  never       death        2     2         671               0.0     28.7198     529.0079              0.0  -0.0071       True
minimum effective-sample-size ratio: 78.8%

What this output tells you. The score-equation check passed: targeting solved the estimating equation of each parameter. The support and nuisance reports are completed, which means that they ran. No row needs attention.

support column

value on this fit

n_followed

671 patients follow never through the second node, the smallest group

max_weight, effective_n

the largest weight is 28.7, for never at the second node, with an effective sample of about 529

share_truncated

0.0 on every row. No probability reached the bound

epsilon

the targeting coefficient. The largest, -0.0679, is on the never-plan readmission row at the second node

converged

targeting converged on every row

The minimum effective-sample-size ratio is 78.8%. These numbers come from the fitted logistic offer and enrollment models. They describe those models, not the population. They do not establish exchangeability, independent censoring, or correct outcome regressions.

Step 10: how much a recorded cause moves the estimate#

A sensitivity analysis asks how strong an unmeasured confounder must be to change the conclusion. No omitted-variable bound is implemented for a longitudinal fit, so this step uses a benchmark. It refits the day-60 differences without one recorded covariate at a time. Each refit uses the same data, method, plans, and protocol. The sensitivity guide describes the approach.

full_design = study.design
benchmarks = []
for dropped in ("age", "baseline_readiness", "identified_needs"):
    reduced = replace(
        full_design,
        baseline=tuple(name for name in full_design.baseline if name != dropped),
        time_varying=tuple(
            tuple(name for name in node if name != dropped) for node in full_design.time_varying
        ),
    )
    refit = (
        CausalStudy(frame, design=reduced, protocol=protocol)
        .identify(levels_by_plan)
        .estimate(method=sequential)
    )
    for cause in ("readmission", "death"):
        names = [incidence("always", cause, 2), incidence("never", cause, 2)]
        moved = refit.contrast(lambda psi: psi[0] - psi[1], names, name=f"{cause} difference")
        full = differences[cause, 2]
        benchmarks.append(
            {
                "dropped": dropped,
                "cause": cause,
                "psi": round(moved.psi, 4),
                "move": round(moved.psi - full.psi, 4),
                "move_in_se": round((moved.psi - full.psi) / full.std_error, 2),
            }
        )
benchmark = pd.DataFrame(benchmarks)
for cause in ("readmission", "death"):
    full = differences[cause, 2]
    print(f"full history, {cause} difference at t=2: {full.psi:.4f}, SE {full.std_error:.4f}")
print(benchmark.to_string(index=False))
full history, readmission difference at t=2: 0.0223, SE 0.0181
full history, death difference at t=2: -0.1376, SE 0.0156
           dropped       cause     psi   move  move_in_se
               age readmission  0.0540 0.0317        1.75
               age       death -0.1276 0.0100        0.64
baseline_readiness readmission  0.0332 0.0109        0.60
baseline_readiness       death -0.1318 0.0057        0.37
  identified_needs readmission  0.0366 0.0143        0.79
  identified_needs       death -0.1201 0.0175        1.12

What this output tells you. Each row drops one covariate. move_in_se divides the move of the day-60 difference by the standard error of the full fit.

dropped covariate

readmission move (SE)

death move (SE)

age

1.75

0.64

baseline_readiness

0.60

0.37

identified_needs

0.79

1.12

Every move is positive here. A dropped covariate leaves both the offer models and the enrollment models, so each move mixes offer confounding with informative censoring. Each covariate sends offers toward patients at higher risk, so its omission makes navigation look worse. The largest death move is 0.0175, against a death difference of -0.1376. The readmission difference is 0.0223, and the move without age is 0.0317. An unrecorded cause that acts like age could therefore change the sign of the readmission difference.

Read the table as a benchmark, not as a bound. An unrecorded cause can be stronger, or act in the other direction.

Method evidence#

The fits on this page use one draw of a synthetic law. The registered studies in the table measure each method over repeated draws. Each study page gives its law, its nuisance models, and the result of each cell.

fit on this page

registered study

law of the study

Step 6, the incidence of each cause

Ordinary competing-risk longitudinal TMLE

another law

Step 7, the event-free difference

Ordinary survival-curve longitudinal TMLE

this page’s law, with both causes merged into one event

Both studies supply the mechanisms exactly. No registered study fits the offer and enrollment models from data, as this page does.

The technical entry lists the evidence and the limits of the method.

What the analysis shows#

The question was whether navigation at discharge and on day 31 lowers readmission by day 60. The estimates describe enrolled adults discharged home, under the plan that offers navigation at both decisions, had every patient stayed enrolled. Each interval is a pointwise 95% interval from the influence curve.

estimand, always minus never at day 60

estimate (95% CI)

reading

readmission incidence

0.0223 (-0.0132, 0.0578)

the data do not show a change. At the 95% level, they rule out a difference below -0.0132 or above 0.0578

death incidence

-0.1376 (-0.1682, -0.1070)

fewer patients die before readmission

event-free survival

0.1153 (0.0748, 0.1558)

more patients reach day 60 alive and not readmitted

Read the readmission result with the death result. A plan that prevents deaths keeps more patients at risk of readmission. A readmission total effect near zero therefore does not show that navigation leaves the readmission hazard unchanged. This analysis does not separate the paths through death. That separation needs a different estimand, with assumptions that this page does not defend.

The analysis does not verify sequential exchangeability, independent censoring, consistency, no interference, or the correctness of the offer and enrollment models. The point estimates rely on those models, because the outcome regressions are misspecified.

the analysis does not support

because

a claim about readmission in a world without death

that is the estimand of Step 8, with more assumptions

a claim about other offer schedules

the fit estimates the always and never plans only

a coverage claim for these intervals

one draw cannot measure coverage

Where to go next#

A program that compares many navigation schedules needs a summary across plans and horizons. MSM projections fit that summary. The examples index lists every tutorial in the program.