MSM projections: summarize three navigation cadences with one trend

A marginal structural model (MSM) projection fits a declared working model to counterfactual means. This notebook fits a straight line in assigned contacts to the means of three navigation cadences. Each step shows its code, its output, and what the output tells you. MSM projections gives the projection, its clever covariate, and the algorithm.

The situation

A hospital network sends adults home after a discharge-home order. At that order, staff assign one of three navigation cadences. Nobody randomized the assignment.

tier	assigned protocol	assigned contacts in 30 days
low	one transition-planning contact	1
medium	planning plus one follow-up contact	2
high	planning plus five follow-up contacts	6

The table separates what the data record from what the question asks about.

quantity	kind	when it is fixed	column
discharge risk, age, baseline support need	observed baseline covariates, W	before assignment	discharge_risk, age, baseline_support_need
assigned cadence	observed treatment, A	at the discharge-home order, which is time zero	cadence
transition score	observed outcome, Y	30 days after discharge	transition_score
the score had the patient been assigned cadence a	counterfactual outcome, Y^a	never for every a. Each patient shows only the score of the assigned cadence	none

Completed contacts, readmission, and death occur after assignment. The data set does not record them, and the analysis must not adjust for them.

The causal graph

discharge_risk ──┬───────────────────► transition_score ◄── baseline_support_need
                 │                       ▲        ▲
                 ▼                       │        │
              cadence ───────────────────┘        │
                 ▲                                │
age ─────────────┴────────────────────────────────┘

arrow	what it asserts	what the analysis does with it
discharge_risk → cadence, age → cadence	staff use risk and age to choose the cadence	with the next row, it makes both variables confounders
discharge_risk → transition_score, age → transition_score	risk and age change the score under every cadence	the adjustment set includes both, which blocks the two back-door paths
baseline_support_need → transition_score	support need changes the score, and it does not change the cadence	it is not a confounder. The adjustment set includes it to explain outcome variance
cadence → transition_score	the effect that the projection summarizes	the only open path from cadence to score after adjustment

absent arrow	what it asserts	what the analysis does with it
no unmeasured cause of both cadence and transition_score	no unmeasured confounding	a scientific judgment. The data cannot check it, and Step 13 asks how strong a violation must be
no arrow from one patient’s cadence to another patient’s score	no interference	a scientific judgment, which the protocol supports with reserved navigator capacity
no node for completed contacts, readmission, or death	these events occur after assignment	completed contacts are a mediator, and death is a competing event. The adjustment set excludes all three. Step 4 states how the protocol handles each

The graph has no censoring node, because every patient has a score. In this synthetic law, the cadence effect does not vary with the covariates, so no variable is an effect modifier. A real program would need to check that claim. The graph says nothing about the shape of the cadence effect. A straight line in contacts is therefore a choice of summary, not a property of the graph.

What you will learn

after this notebook you can	the step that shows it
say why unadjusted cadence means mislead here	Step 3
estimate one counterfactual mean per cadence	Step 6
write the board’s question as a projection and fit it	Step 7
code the cadences for the question the board asks	Step 8
test whether the line fits the three means	the failure mode, Step 9
say why the projection weight is part of the estimand	Step 10
check the projection against the per-cadence report	Step 11
read the support, truncation, and sensitivity evidence	Steps 12 and 13

The causal question and the estimand

The program board wants one number for its report: how the mean 30-day score changes with assigned contacts across the three cadences.

element	definition
target population	adults with a discharge-home order at a participating hospital. Step 4 records the eligibility
interventions	assign every eligible patient low, medium, or high, analyzed as assigned
counterfactual means	E[Y^a] for each cadence a
estimand	the coefficients (β0, β1) of the line β0 + β1 c(a) that minimizes the sum over cadences of h(a) (E[Y^a] − β0 − β1 c(a))². Here c(a) is the assigned contacts and h(a) = 1

The slope β1 has a value whether or not the three means lie on a line. When they do not, β1 is a summary of the means. It is not the effect of one more contact.

simpler estimand	why it does not answer the board’s question
the difference in observed mean scores	confounding biases it. Step 3 shows the observed order of two cadences reversed
the three means, or the two contrasts against low	they answer one question per cadence and give no single per-contact number. They stay the primary result, and Step 11 reports them
the contact coefficient of an outcome regression	it is a projection of the observed data, so its value depends on the covariate model. It does not target E[Y^a]
the slope of a regression weighted by stabilized inverse probabilities	it uses h(a) = P(A = a). That is a different estimand when the cadence shares differ and the line does not fit (Step 10)

Identification needs four assumptions.

assumption	what it means for this program	basis
consistency	each cadence always means its declared script and contacts	scientific judgment
no interference	one patient’s cadence does not change another patient’s score	scientific judgment
no unmeasured confounding	risk, age, and support need block every common cause of cadence and score	scientific judgment
positivity	each covariate profile has some chance of each cadence	checkable in part, through the support report of Step 12

Comparisons on this page

Each comparison isolates one reason why two numbers differ. The page names that reason before the output.

step	comparison	estimands	mechanism that makes them differ	evidence
3	observed cadence means against the true means	observed conditional means against E[Y^a]	risk and age confound the assignment	the confounder means by cadence
8	the contact coding against the step coding	two projections with different designs	each coding spaces the cadences differently, and the line does not fit	the population slope of each coding
9	the fitted line against the three cadence means	projection values against E[Y^a]	the gain per contact falls from low to high	the miss at medium, with its interval
10	the uniform weight against the share and 1:10:1 weights	three projections with different weights	the line does not fit, so the weight decides which cadences it fits best	the population slope under each weight
11	the saturated projection against the cadence contrasts	the same estimand, written two ways	none. The two must agree	equal estimates and influence curves

Why this method

Each method choice answers one feature of the problem.

feature of the problem	method choice	reason
one assignment at time zero, with baseline confounders only	point-treatment TMLE	no later treatment decision, time-varying confounder, or censoring needs longitudinal TMLE
a treatment with three unordered labels	a multinomial treatment model g	each cadence gets its own probability, and positivity applies to each cadence
a board question about one per-contact summary	an MSM projection with a known design and a uniform weight	it reduces three means to two coefficients, each with an influence curve and an interval
a summary whose working model can be wrong	a lack-of-fit contrast and a saturated control	Steps 9 and 11 show what the line hides
a score observed for every patient	no missing-outcome model	no score is missing in this data set
finite-dimensional learners and no documented score range	nuisances fitted in sample, without cross-fitting	Step 6 gives the reason

The projection has three costs. Its coefficients mean what the working model says. You code each cadence label as a number yourself. The slope is not the effect of one more contact unless the three means lie on a line.

term	plain meaning
estimand	the number the question asks for, written before any model is chosen. See estimands
counterfactual mean	the mean score if every eligible patient were assigned one cadence
working model	the line you declare, as a design with named terms. It summarizes the means and does not claim to be the true response
projection weight	how much each cadence counts when the line is fitted to the means. Here every cadence counts equally
nuisance	a model the estimate needs but the question does not ask about. Here, the outcome regression Q and the treatment mechanism g. See point-treatment TMLE
saturated model	a working model with one free coefficient per cadence, so it reproduces the three means exactly

Step 1: set up

The setup imports the learners and prints the installed cleverly version. Every fit below passes linear learners, a fold count, and a random seed. A rerun therefore reproduces the stored outputs.

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

import cleverly

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

cleverly 0.1.3

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

Step 2: the data

make_multi_arm draws a three-arm confounded law with known counterfactual means. The code renames its columns to the program’s names. It records the contact count of each cadence, prints the first rows, and prints the known means.

from cleverly.datasets import make_multi_arm

frame, truth = make_multi_arm(n=3_000, seed=61)
frame = frame.rename(
    columns={
        "Y": "transition_score",
        "A": "cadence",
        "W1": "discharge_risk",
        "W2": "age",
        "W3": "baseline_support_need",
    }
)
ARMS = ("low", "medium", "high")
CONTACTS_30D = {"low": 1.0, "medium": 2.0, "high": 6.0}
contacts = np.array([CONTACTS_30D[arm] for arm in ARMS])
population = np.array([truth[f"ey[{arm}]"] for arm in ARMS])

print("rows and columns:", frame.shape)
print("cadence labels:", sorted(frame["cadence"].unique()))
print(frame.head().round(3))
print()
print("known counterfactual means of the synthetic law:")
for arm, value in zip(ARMS, population, strict=True):
    print(f"  ey[{arm}] ({CONTACTS_30D[arm]:.0f} contacts): {value:.4f}")
gain_per_contact = np.diff(population) / np.diff(contacts)
print(f"gain per contact, low to medium:  {gain_per_contact[0]:.2f}")
print(f"gain per contact, medium to high: {gain_per_contact[1]:.2f}")

rows and columns: (3000, 5)
cadence labels: ['high', 'low', 'medium']
   transition_score cadence  discharge_risk    age  baseline_support_need
0             1.266    high          -0.567 -1.057                  1.062
1            -0.691     low          -1.129 -1.779                  0.655
2             1.325  medium           0.841 -0.949                 -0.565
3             5.003  medium           2.780 -1.017                 -0.506
4             0.847  medium           1.179  0.745                  1.612

known counterfactual means of the synthetic law:
  ey[low] (1 contacts): -0.0000
  ey[medium] (2 contacts): 0.6000
  ey[high] (6 contacts): 1.4400
gain per contact, low to medium:  0.60
gain per contact, medium to high: 0.21

What this output tells you. Each row is one discharge, and cadence holds a text label. The published truths use the generator’s labels, so the labels stay unchanged.

The generator draws the three covariates as independent standard-normal variables. Risk and age enter both the cadence probabilities and the score. Support need enters only the score. These are the arrows of the causal graph. The score has no unit and no documented range.

The true means are 0 for low (printed as -0.0000 after rounding), 0.6000 for medium, and 1.4400 for high. The gain per contact falls from 0.60 between low and medium to 0.21 between medium and high. No straight line in assigned contacts passes through the three means, so the working model of Step 7 is misspecified on purpose.

A real program has no truth. Every comparison against it below is a teaching device.

Step 3: association first

A confounder is a variable that changes both the assigned cadence and the outcome. The code compares the cadences before any adjustment. It prints the observed mean score, the means of the two confounders, the share of patients, and the true counterfactual mean.

observed = frame.groupby("cadence")[["transition_score", "discharge_risk", "age"]].mean()
observed = observed.loc[list(ARMS)]
observed["share"] = frame["cadence"].value_counts(normalize=True).loc[list(ARMS)]
observed["true mean"] = population
print(observed.round(3))

         transition_score  discharge_risk    age  share  true mean
cadence                                                           
low                -0.082          -0.132 -0.128  0.272      -0.00
medium              1.324           0.532 -0.343  0.368       0.60
high                0.808          -0.387  0.474  0.360       1.44

What this output tells you. The observed means put medium above high, at 1.324 against 0.808. The true means reverse that order, at 0.60 against 1.44.

The confounded assignment causes the reversal. The table compares the two confounders in the medium and high groups.

confounder	medium mean	high mean	effect on the score in this law
discharge_risk	0.532	-0.387	a higher value raises the score
age	-0.343	0.474	a higher value lowers the score

Both confounders push the observed medium mean above the observed high mean. A trend fitted to these observed means would summarize the confounding too.

The shares are 0.272, 0.368, and 0.360. Step 10 uses them to show that the projection weight is part of the estimand.

Step 4: write the protocol

A StudyProtocol records the scientific design before any model runs. The record gets a fingerprint, and every result fitted from it carries that fingerprint. This page starts from navigation_protocol(), the protocol of the shared study design. dataclasses.replace lists the three cadences as three strategies, each with its version.

from dataclasses import replace

from cleverly.datasets import navigation_protocol

program = navigation_protocol()
protocol = replace(
    program,
    outcome="Patient-reported transition score in synthetic units, with no documented range",
    time_zero="Discharge-home order, after baseline measurement and before cadence assignment",
    treatment_strategies=(
        "Assign the low navigation cadence",
        "Assign the medium navigation cadence",
        "Assign the high navigation cadence",
    ),
    treatment_versions=(
        "One transition-planning contact within 30 days",
        "Transition planning plus one follow-up contact within 30 days",
        "Transition planning plus five follow-up contacts within 30 days",
    ),
    intercurrent_event_handling=(
        program.intercurrent_event_handling[0],
        "Analyze the assigned cadence regardless of completed contacts",
        "The score has no worst value, so the protocol handles death before day 30 "
        "by the hypothetical strategy: the score had the patient not died",
    ),
    assumption_rationale=(
        "The recorded baseline variables cover the measured common causes of cadence and score",
        "Every cadence uses the same script and access rules, which supports consistency",
        program.assumption_rationale[2],
    ),
)
print("\n".join(protocol.summary_lines()))

causal study protocol: schema 1; ec77f97f713c2306
target population: Adults with a discharge-home order at a participating hospital during the enrollment period
eligibility: ['Age 18 years or older', 'Discharge home ordered at a participating hospital']
time zero: Discharge-home order, after baseline measurement and before cadence assignment
treatment strategies: ['Assign the low navigation cadence', 'Assign the medium navigation cadence', 'Assign the high navigation cadence']
treatment versions: ['One transition-planning contact within 30 days', 'Transition planning plus one follow-up contact within 30 days', 'Transition planning plus five follow-up contacts within 30 days']
outcome: Patient-reported transition score in synthetic units, with no documented range
horizon: 30 days after discharge
intercurrent-event handling: ['Use the transition score regardless of readmission', 'Analyze the assigned cadence regardless of completed contacts', 'The score has no worst value, so the protocol handles death before day 30 by the hypothetical strategy: the score had the patient not died']
interference unit: Individual patient
assumption rationale: ['The recorded baseline variables cover the measured common causes of cadence and score', 'Every cadence uses the same script and access rules, which supports consistency', 'Reserved navigator capacity and access controls support no interference']

What this output tells you. The first line gives the schema version and the fingerprint ec77f97f713c2306. The other lines repeat each field. This page changes six fields of navigation_protocol(), and point-treatment TMLE reads the others.

protocol field	what this page writes
outcome	the score is in synthetic units and has no documented range. The program records a share of the maximum score
time zero	the discharge-home order, before cadence assignment
treatment strategies and versions	the three cadences, with their contacts in words
intercurrent-event handling	the assigned cadence is analyzed regardless of completed contacts. The program scores death as the worst score, and this score has no worst value. This page uses the hypothetical strategy for death
assumption rationale	the common causes of cadence and score, and one script for every cadence

Under the hypothetical strategy, a real analysis treats death as censoring. It then needs exchangeability and positivity for death given the recorded history. This synthetic law draws no deaths, so the choice does not change the estimate here.

A score with no documented range has no range to declare as q_bounds. Step 6 uses that fact.

The protocol has no field for the contact mapping, the working model, or the projection weight. The typed estimand in Step 7 owns all three.

Step 5: design and identification

The design is an ordinary point-treatment design with a three-level treatment. CounterfactualMean asks for one mean per cadence. The projection in Step 7 summarizes those means.

from cleverly import CausalStudy, CounterfactualMean, PointTreatment

study = CausalStudy(
    frame,
    design=PointTreatment(
        outcome="transition_score",
        treatment="cadence",
        adjustment=("discharge_risk", "age", "baseline_support_need"),
    ),
    protocol=protocol,
)
arm_means_effect = study.identify(CounterfactualMean())
print(arm_means_effect.summary(protocol="fingerprint"))

counterfactual mean under each treatment, E[Y^a]
identified by explicit-adjustment: E_W[E(transition_score | cadence=a, W)] for a in ['high', 'low', 'medium']
adjustment/history: ['discharge_risk', 'age', 'baseline_support_need']
required nuisances: ['outcome_regression', 'treatment_mechanism']
assumptions:
  - consistency: Y = Y^a when A = a
  - no interference: one unit's potential outcome does not depend on other units' treatment assignments
  - no unmeasured confounding: Y^a is independent of A given W
  - positivity: P(cadence = a | W) > 0 almost surely for every supported treatment level a in ['high', 'low', 'medium']
causal study protocol: schema 1; ec77f97f713c2306

What this output tells you. The formula averages the outcome regression over all patients, once for each cadence. The required nuisances are the outcome regression and the treatment mechanism. The four assumptions are the ones that the causal question lists with their basis.

The positivity condition names high, low, and medium. Each level needs a positive probability at every covariate value where the target places mass. The summary then prints the fingerprint of the stored protocol.

Step 6: estimate each cadence’s mean

Start with the per-cadence report, because the projection is a summary of it. The linear outcome model and the multinomial logistic treatment model are correctly specified for this law.

setting	value	what it does
outcome_learner	LinearRegression()	fits Q, the expected score given the cadence and the covariates
treatment_learner	LogisticRegression(max_iter=1000, random_state=61)	fits g, the probability of each cadence given the covariates
CrossFitting(enabled=False)	no outer folds	fits each nuisance on every row
Runtime(random_state=61, n_jobs=1)	fixed seed, one process	makes the fit reproducible

The fit runs in sample for two reasons. A cross-fitted fit of a continuous outcome declares the outcome’s range as q_bounds, and this score has no documented range. The learners are finite-dimensional regressions. That makes plausible the data-reuse condition that CV-TMLE removes. The cross-fitting tutorial shows the fold layer on a score with a recorded range.

Targeting makes a small update to Q, weighted by g, that removes first-order bias. Targeting and bounds defines it.

from cleverly import CrossFitting, ModelSpec, Runtime, TMLEMethod

method = TMLEMethod(
    models=ModelSpec(
        outcome_learner=LinearRegression(),
        treatment_learner=LogisticRegression(max_iter=1000, random_state=61),
    ),
    cross_fitting=CrossFitting(enabled=False),
    runtime=Runtime(random_state=61, n_jobs=1),
)
arm_result = arm_means_effect.estimate(method=method)

arm_table = arm_result.to_frame().set_index("estimand").loc[[f"ey[{arm}]" for arm in ARMS]]
arm_table = arm_table[["psi", "std_err", "ci_lower", "ci_upper"]]
arm_table["true mean"] = population
print(arm_table.round(4))
print()
print("protocol fingerprint on the fit:", arm_result.provenance.protocol_fingerprint)

               psi  std_err  ci_lower  ci_upper  true mean
estimand                                                  
ey[low]    -0.0038   0.0365   -0.0754    0.0678      -0.00
ey[medium]  0.6399   0.0386    0.5641    0.7156       0.60
ey[high]    1.4642   0.0390    1.3877    1.5408       1.44

protocol fingerprint on the fit: ec77f97f713c2306

What this output tells you. Each row is one counterfactual mean. The medium estimate is 0.6399 with the interval (0.5641, 0.7156), and the true mean is 0.60. The high estimate is 1.4642, and the low estimate is -0.0038.

Each interval contains its true mean on this draw. That is one draw, not a coverage result. The interval comes from the influence curve, which measures how much each row moves the estimate. See inference. The fit carries the protocol fingerprint ec77f97f713c2306.

Step 7: declare the working model and fit the trend

The design function receives one cadence label and the covariate frame. It returns the design matrix for that cadence, with one column per term. The terms tuple names the columns.

The declaration leaves weights unset, which sets the projection weight h(a, V) to 1 for every cadence. The same method fits the projection.

from cleverly import MSMProjection
from cleverly.msm import MSM


def contacts_design(arm, data):
    """One row per patient: an intercept and the cadence's assigned contact count."""
    return np.column_stack([np.ones(len(data)), np.full(len(data), CONTACTS_30D[arm])])


trend = MSM(design=contacts_design, terms=("(intercept)", "assigned contacts"), design_kind="known")
trend_effect = study.identify(MSMProjection(trend))
trend_result = trend_effect.estimate(method=method)
print(trend_result.summary(protocol="fingerprint"))

Targeted maximum likelihood estimation
======================================
n = 3000; covariates = 3; arm shares: high=0.36, low=0.272, medium=0.368
causal estimand: projection of counterfactual means onto a working model
identification: explicit-adjustment; identified msm-indexed plug-in functional of E(transition_score | cadence=a, W) for a in ['high', 'low', 'medium'] with influence correction
required nuisances: outcome_regression, treatment_mechanism
identification assumptions: consistency: Y = Y^a when A = a; no interference: one unit's potential outcome does not depend on other units' treatment assignments; no unmeasured confounding: Y^a is independent of A given W; positivity: P(cadence = a | W) > 0 almost surely for every supported treatment level a in ['high', 'low', 'medium']; the working model and its weights are known functions of (a, V): neither phi nor h depends on the observed-data law, so the influence function carries no term for estimating them; the weighted Gram matrix is invertible, so the projection is a single coefficient vector rather than a set of them
causal study protocol: schema 1; ec77f97f713c2306
outcome family: gaussian; fluctuation: logistic (iterative)
TMLE (in-sample nuisances): nuisances fitted in-sample (cross_fit=False)
propensity truncated to [0.0114, 0.9886]; ATT/ATC to [0.025, 0.975]
outcome scaled from [-6.044, 7.723] to [0, 1]
data 248bf6e493a8a2be | folds 9ba86d4d13cb065e | cleverly 0.1.3 | 2026-10-16T01:17:23+00:00

estimand                psi      std_err   95% CI                p_value  support     
----------------------  -------  --------  --------------------  -------  ------------
msm[(intercept)]        -0.1068  0.03703   [-0.17937, -0.03423]  0.0039   experimental
msm[assigned contacts]  0.26875  0.008594  [0.25191, 0.2856]     <1e-4    experimental

experimental: No registered study covers this configuration. Nearest row: T18 (point-msm), differing on treatment.

simultaneous 95% bands (multiplier bootstrap, critical value 2.107 vs 1.960 pointwise) [experimental]:
  msm[(intercept)] [-0.18484, -0.028769]
  msm[assigned contacts] [0.25064, 0.28686]

What this output tells you. The identification assumptions line lists the four causal assumptions from Step 5 and adds two.

added assumption	what it means here
the working model and its weights are known functions	contacts_design and the uniform weight do not depend on the data. design_kind="known" declares this for contacts_design. The uniform weight needs no declaration
the weighted Gram matrix is invertible	the three cadences have distinct contact counts, so one line is closest

The configuration lines have the same meaning as in point-treatment TMLE. Here the fluctuation solves one score equation per coefficient, not one per cadence.

The slope is 0.26875 transition-score units per additional assigned contact. The intercept is -0.1068, the line’s value at zero contacts, which no cadence assigns.

The simultaneous 95% bands are built to contain both coefficients at the same time, so their critical value is 2.107 instead of 1.960. Quote the pointwise interval when you report the slope alone.

The support column marks both coefficients and the bands experimental, because no registered study covers this configuration. The support label makes no coverage claim. The review probe reviews/notebook-review/probes/msm-projections-final/sweep.py repeats this fit on 4,000 draws of this law at n = 3,000. The interval covered on 3,762 of 4,000 draws, a rate of 0.9405. The probe does not isolate the cause of the shortfall.

With a uniform weight, the slope is a fixed contrast of the three means: (-2 E[Y(low)] - E[Y(medium)] + 3 E[Y(high)]) / 14. Step 9 checks that identity on the population means.

Step 8: code the cadences for the question the board asks

The contact counts are the program’s decision, not the estimator’s. The number that a design function writes for each label sets the question the slope answers. The code fits a second known design that codes each cadence by its step, 0, 1, and 2. It prints both slopes beside the population slope of each coding.

CADENCE_STEP = {"low": 0.0, "medium": 1.0, "high": 2.0}


def cadence_step_design(arm, data):
    """One row per patient: an intercept and the cadence's step, 0, 1, or 2."""
    return np.column_stack([np.ones(len(data)), np.full(len(data), CADENCE_STEP[arm])])


step_trend = MSM(
    design=cadence_step_design, terms=("(intercept)", "cadence step"), design_kind="known"
)
step_result = study.identify(MSMProjection(step_trend)).estimate(method=method)

codings = (
    ("assigned contacts", CONTACTS_30D, trend_result),
    ("cadence step", CADENCE_STEP, step_result),
)
coding_rows = []
for term, coding, result in codings:
    values = np.array([coding[arm] for arm in ARMS])
    coding_line, *_ = np.linalg.lstsq(np.column_stack([np.ones(3), values]), population, rcond=None)
    slope = result[f"msm[{term}]"]
    coding_rows.append(
        {
            "coefficient": f"msm[{term}]",
            "low, medium, high": ", ".join(f"{value:.0f}" for value in values),
            "slope": f"{slope.psi:.4f}",
            "std_err": f"{slope.std_error:.4f}",
            "95% CI": f"({slope.ci[0]:.4f}, {slope.ci[1]:.4f})",
            "population slope": f"{coding_line[1]:.4f}",
        }
    )
print(pd.DataFrame(coding_rows).to_string(index=False))

           coefficient low, medium, high  slope std_err           95% CI population slope
msm[assigned contacts]           1, 2, 6 0.2688  0.0086 (0.2519, 0.2856)           0.2657
     msm[cadence step]           0, 1, 2 0.7341  0.0224 (0.6902, 0.7779)           0.7200

What this output tells you. The two rows read the same three means through two codings.

coefficient	slope	interval	population slope	question it answers
msm[assigned contacts]	0.2688	(0.2519, 0.2856)	0.2657	the change per assigned contact
msm[cadence step]	0.7341	(0.6902, 0.7779)	0.7200	the change per cadence step

The step coding counts the move from medium to high as one unit. The contact coding counts the same move as four contacts. Neither line passes through the three means. The board asks for the change per contact, so the rest of this page reads the contact slope. Each interval contains its population slope on this draw.

Under the uniform weight, the per-step slope is half the difference between the high and low means. Its population value is 1.44 / 2 = 0.7200. Step 11 prints the fitted high vs low contrast as 1.4681, and half of that is 0.734.

Write the coding in the design function, as both fits here do. The labels have no order of their own. Sorted as text, they are high, low, and medium.

Step 9: the failure mode, a projection is not a fitted curve

The working model from Step 7 is wrong for this law. This step compares two estimands: the projection’s value at each cadence, β0 + β1 c(a), and the cadence mean E[Y^a]. The mechanism that separates them is the falling gain per contact that Step 2 prints. The code compares the line with the means, first in the population and then in the fits. It measures the miss at medium with arm_result.contrast(...), which reads the joint influence curve of the three cadence means. It also computes the population slope as the fixed contrast of the three means.

design = np.vstack([contacts_design(arm, frame.iloc[:1]) for arm in ARMS])
names = [f"msm[{term}]" for term in trend.terms]
projection, *_ = np.linalg.lstsq(design, population, rcond=None)
estimated_line = design @ np.array([trend_result[name].psi for name in names])

in_population = pd.DataFrame(
    {
        "contacts": contacts,
        "population mean": population,
        "population line": design @ projection,
        "line minus mean": design @ projection - population,
    },
    index=list(ARMS),
)
in_fits = (
    arm_table[["psi", "ci_lower", "ci_upper"]]
    .set_axis(list(ARMS))
    .rename(columns={"psi": "estimated mean"})
)
in_fits["estimated line"] = estimated_line
in_fits["line minus mean"] = estimated_line - in_fits["estimated mean"]
print("in the population:")
print(in_population.round(4).to_string())
print()
print("in the fits:")
print(in_fits.round(4).to_string())
step7_medium_miss = in_fits.loc["medium", "line minus mean"]
print(f"absolute Step 7 miss at medium: {abs(step7_medium_miss):.4f}")
print()

# The least-squares line at medium is a fixed weighted sum of the three means.
misfit_weights = (design @ np.linalg.pinv(design))[1] - np.eye(3)[1]
medium_misfit = arm_result.contrast(
    lambda means: float(misfit_weights @ means),
    [f"ey[{arm}]" for arm in ARMS],
    name="medium: line minus mean",
    gradient=lambda means: misfit_weights,
)
low, high = medium_misfit.ci
print(f"medium: line minus mean, {medium_misfit.psi:.4f}  CI=({low:.4f}, {high:.4f})")
print(f"absolute difference from the Step 7 miss: {abs(medium_misfit.psi - step7_medium_miss):.4f}")
print()
print(f"population projection: intercept {projection[0]:.4f}, slope {projection[1]:.4f}")
contrast_slope = (-2 * population[0] - population[1] + 3 * population[2]) / 14
print(f"slope as the fixed contrast of the three means: {contrast_slope:.4f}")
for name, target in zip(names, projection, strict=True):
    low, high = trend_result[name].ci
    print(f"{name}: CI=({low:.4f}, {high:.4f})  population projection {target:.4f}")

in the population:
        contacts  population mean  population line  line minus mean
low          1.0            -0.00           0.1486           0.1486
medium       2.0             0.60           0.4143          -0.1857
high         6.0             1.44           1.4771           0.0371

in the fits:
        estimated mean  ci_lower  ci_upper  estimated line  line minus mean
low            -0.0038   -0.0754    0.0678          0.1620           0.1658
medium          0.6399    0.5641    0.7156          0.4307          -0.2092
high            1.4642    1.3877    1.5408          1.5057           0.0415
absolute Step 7 miss at medium: 0.2092

medium: line minus mean, -0.2084  CI=(-0.2566, -0.1602)
absolute difference from the Step 7 miss: 0.0008

population projection: intercept -0.1171, slope 0.2657
slope as the fixed contrast of the three means: 0.2657
msm[(intercept)]: CI=(-0.1794, -0.0342)  population projection -0.1171
msm[assigned contacts]: CI=(0.2519, 0.2856)  population projection 0.2657

What this output tells you. In the population, the line misses medium by about 0.19, since line minus mean is -0.1857. It overshoots low by 0.1486 and high by 0.0371.

In the fits, the absolute Step 7 miss at medium line reports 0.2092. The medium: line minus mean row measures that miss with an interval.

part of the row	value on this draw	what it means
estimate	-0.2084	the line fitted to the three estimated means, minus the estimated medium mean
interval	(-0.2566, -0.1602)	it excludes zero, so the data reject, at the 5% level, a line through all three means
population miss	-0.1857	the interval contains it

With three cadences and two coefficients, a nonzero miss at one cadence is equivalent to a nonzero miss at every cadence. The population misses stand in the ratio 4 : -5 : 1. The contrast call uses the joint influence curve of the three means, so their correlation enters the interval. Its estimate differs from the Step 7 miss by 0.0008 on this draw, because the cadence fit and the MSM fit target separately.

check	value on this draw
population projection	intercept -0.1171, slope 0.2657
slope as the fixed contrast of the means	0.2657
slope interval from the fit	(0.2519, 0.2856), which contains 0.2657
intercept interval from the fit	(-0.1794, -0.0342), which contains -0.1171

Read the slope as “the best linear summary of the cadence response under a uniform weight”. Do not read it as “the causal effect of one more contact”. No patient was assigned more than six contacts, so the projection says nothing about a cadence with ten.

Step 10: the weight is part of the estimand

The projection weight sets how much each cadence counts when the line is fitted. The code fits the population line again, with each cadence weighted by its share from Step 3. It then declares a fixed weight of 1, 10, and 1 for low, medium, and high. It prints the population slope under that weight, and it fits the projection with the weight declared as weights_kind="known".

shares = observed["share"].to_numpy()
# Weighted least squares: solve (X' W X) b = X' W y, with W the diagonal matrix of weights.
share_projection = np.linalg.solve(
    design.T @ (shares[:, None] * design), design.T @ (shares * population)
)
slope = trend_result["msm[assigned contacts]"]
print(
    "cadence shares:",
    ", ".join(f"{arm} {share:.3f}" for arm, share in zip(ARMS, shares, strict=True)),
)
print(f"population slope, uniform weight: {projection[1]:.4f}")
print(f"population slope, share weight:   {share_projection[1]:.4f}")
print(f"shift from the share weight:      {abs(share_projection[1] - projection[1]):.4f}")
print(f"standard error of the Step 7 slope: {slope.std_error:.4f}")
print()

# A fixed weight that puts most of the projection on medium.
FIXED_WEIGHT = {"low": 1.0, "medium": 10.0, "high": 1.0}
fixed = np.array([FIXED_WEIGHT[arm] for arm in ARMS])
fixed_projection = np.linalg.solve(
    design.T @ (fixed[:, None] * design), design.T @ (fixed * population)
)
weighted_trend = MSM(
    design=contacts_design,
    terms=trend.terms,
    weights=lambda arm, data: np.full(len(data), FIXED_WEIGHT[arm]),
    weights_kind="known",
    design_kind="known",
)
weighted_result = study.identify(MSMProjection(weighted_trend)).estimate(method=method)
weighted_slope = weighted_result["msm[assigned contacts]"]
print(f"population slope, 1:10:1 weight:  {fixed_projection[1]:.4f}")
for label, estimate in (("uniform", slope), ("1:10:1", weighted_slope)):
    low, high = estimate.ci
    print(f"fitted slope, {label} weight: {estimate.psi:.4f}  CI=({low:.4f}, {high:.4f})")

cadence shares: low 0.272, medium 0.368, high 0.360
population slope, uniform weight: 0.2657
population slope, share weight:   0.2593
shift from the share weight:      0.0065
standard error of the Step 7 slope: 0.0086

population slope, 1:10:1 weight:  0.2400
fitted slope, uniform weight: 0.2688  CI=(0.2519, 0.2856)
fitted slope, 1:10:1 weight: 0.2398  CI=(0.2214, 0.2582)

What this output tells you. The uniform weight gives a population slope of 0.2657. The share weight gives 0.2593. The means are the same, so only the weight moved the slope.

The shift is 0.0065, and the standard error of the Step 7 slope is 0.0086. That standard error describes one slope, not the difference between two slopes. This illustration uses population means and sample shares. It does not test whether two fitted slopes differ.

The 1:10:1 weight moves the population slope to 0.2400. The fitted slope under that weight is 0.2398, with the interval (0.2214, 0.2582). On this draw, each interval contains its own population slope and not the other.

fitted slope	interval	contains 0.2657	contains 0.2400
uniform weight, 0.2688	(0.2519, 0.2856)	yes	no
1:10:1 weight, 0.2398	(0.2214, 0.2582)	no	yes

A change of weight or terms changes the estimand, not only the estimate. A misspecified working model still defines a parameter. Interval validity still needs the causal and nuisance conditions.

Declare a weights= function as known only when it does not depend on the data. The 1:10:1 weight qualifies. The sample shares and a function of the fitted treatment mechanism do not. Such a weight, declared as known, still fits. That fit omits the weight’s influence-curve term, so its standard error can be too small or too large. The technical entry names the term that an estimated weight or design adds.

Step 11: the control, a saturated working model

A saturated model has one free coefficient per cadence. It cannot be misspecified, so it represents the three means exactly. The code fits it beside the average treatment effect (ATE) contrasts against low, with the same method. Each row pairs one coefficient with its counterpart.

from cleverly import ATE


def cadence_indicators(arm, data):
    """One free coefficient per cadence: an intercept and two cadence indicators."""
    n = len(data)
    return np.column_stack(
        [
            np.ones(n),
            np.full(n, float(arm == "medium")),
            np.full(n, float(arm == "high")),
        ]
    )


saturated = MSM(
    design=cadence_indicators,
    terms=("(intercept)", "medium vs low", "high vs low"),
    design_kind="known",
)
saturated_result = study.identify(MSMProjection(saturated)).estimate(method=method)
arm_contrasts = study.identify(ATE(reference="low")).estimate(method=method)

pairs = (
    ("msm[(intercept)]", arm_result, "ey[low]"),
    ("msm[medium vs low]", arm_contrasts, "ate[medium vs low]"),
    ("msm[high vs low]", arm_contrasts, "ate[high vs low]"),
)
rows = []
for name, counterpart, key in pairs:
    coefficient, match = saturated_result[name], counterpart[key]
    rows.append(
        {
            "coefficient": name,
            "psi": coefficient.psi,
            "ci": f"({coefficient.ci[0]:.4f}, {coefficient.ci[1]:.4f})",
            "counterpart": key,
            "counterpart psi": match.psi,
            "counterpart ci": f"({match.ci[0]:.4f}, {match.ci[1]:.4f})",
            "influence curves agree within 1e-12": bool(
                np.allclose(coefficient.influence_curve, match.influence_curve, rtol=0, atol=1e-12)
            ),
        }
    )
print(pd.DataFrame(rows).round(4).to_string(index=False))

       coefficient     psi                ci        counterpart  counterpart psi    counterpart ci  influence curves agree within 1e-12
  msm[(intercept)] -0.0038 (-0.0754, 0.0678)            ey[low]          -0.0038 (-0.0754, 0.0678)                                 True
msm[medium vs low]  0.6437  (0.5562, 0.7312) ate[medium vs low]           0.6437  (0.5562, 0.7312)                                 True
  msm[high vs low]  1.4681  (1.3805, 1.5557)   ate[high vs low]           1.4681  (1.3805, 1.5557)                                 True

What this output tells you. The intercept is -0.0038, the same value as the low mean in Step 6. The other coefficients are 0.6437 and 1.4681, the same as the medium versus low and high versus low contrasts. Each pair of influence curves agrees within 1e-12 at every row, so the intervals match too.

A saturated projection is a reparameterization of the per-cadence report, not a different analysis. It is also the fallback for a board that does not want to commit to a contact spacing.

A saturated check cannot see a component that vanishes when the model fits. The projection weight and a link’s curvature term are two examples. The technical entry describes the checks that do test them.

Step 12: diagnostics, what the fit can check

include_retargets=True runs the truncation curve with the assessment. Truncation clips each fitted cadence probability into a bound, so no row gets an extreme weight. The curve retargets each coefficient at a range of bounds without refitting the nuisance models. The support report describes positivity in the fitted treatment model. See diagnostics.

assessment = trend_result.assess(include_retargets=True)
ledger = assessment.to_frame()
ran = ledger[~ledger["status"].isin(["unavailable", "not_applicable"])]
print(ran[["surface", "check", "status"]].to_string(index=False))
print("needs attention:", tuple(item.name for item in assessment.attention))

support = assessment.report("support")
scores = assessment.report("score_equations")
curve = assessment.report("truncation_curve")
print()
print(support.summary())
print()
print(scores.summary())
print()
slope_curve = curve.loc[
    curve["estimand"] == "msm[assigned contacts]",
    ["bound", "psi", "delta_from_fitted", "ci_lower", "ci_upper", "truncated_fraction"],
]
print(slope_curve.round(4).to_string(index=False))
print(f"largest absolute slope movement: {slope_curve['delta_from_fitted'].abs().max():.4f}")

    surface            check    status
 validation  score_equations    passed
 validation          support   warning
 validation  nuisance_models completed
diagnostics truncation_curve completed
diagnostics           refute  deferred
needs attention: ('support',)

Positivity / overlap diagnostics
--------------------------------
n = 3000; propensity truncated to [0.0114, 0.9886]

g(W) quantile  0%      1%      5%      25%     50%     75%     95%     99%     100%  
-------------  ------  ------  ------  ------  ------  ------  ------  ------  ------
g[high]        0.0006  0.0153  0.0478  0.1684  0.3291  0.5263  0.7641  0.8681  0.9546
g[low]         0.0148  0.0799  0.1332  0.2191  0.2755  0.3260  0.4022  0.4494  0.5720
g[medium]      0.0044  0.0168  0.0470  0.1703  0.3416  0.5430  0.7773  0.8897  0.9847

arm     n     effective n  ESS / n  top 1% weight  top 5% weight
------  ----  -----------  -------  -------------  -------------
high    1079  514.9        0.477    0.086          0.218        
low     816   727.2        0.891    0.032          0.105        
medium  1105  655.3        0.593    0.068          0.190        

threshold  P(g < t)  P(g > 1-t)
---------  --------  ----------
0.01       0.0033    0.0000    
0.025      0.0136    0.0001    
0.05       0.0371    0.0009    
0.1        0.0963    0.0041    

truncated: 36 unit(s) (1.20%); most extreme untruncated g(W) = 0.011537

group  equation  target rows  Kish-equivalent rows  Kish / target  Kish / total  top 1% load  top 5% load  max |w h|  zero load  g bound           clipped   
-----  --------  -----------  --------------------  -------------  ------------  -----------  -----------  ---------  ---------  ----------------  ----------
msm    h_msm1    3000/3000    1032.8                0.344          0.344         0.105        0.261        268.5      0          [0.0114, 0.9886]  36 (1.20%)
(absolute residual-multiplier load from the fitted equation; Kish values describe residual-multiplier concentration, not residual contributions or estimator information, influence, or precision; target rows are score-mask rows and include structural zeros)
max |clever covariate| (msm): 268.5

VERDICT: some truncation. 1.2% of units were clipped at the marginal treatment-mechanism bound [0.0114, 0.9886]. Report truncation_curve() alongside the estimate so readers can see how much the answer depends on it. The narrowest arm has a Kish-equivalent weight count of 48% of its nominal rows. This is concentration, not estimator information, and no threshold is applied because none is derived. Absolute-load concentration is greatest for group 'msm', equation 'h_msm1': 1032.8 Kish-equivalent rows out of 3000 score-mask rows (34%; 34% of 3000 total rows). This is a descriptive residual-multiplier concentration index. It excludes the residual and is not estimator effective sample size, information, actual influence, precision, or a positivity test.

Score-equation check
--------------------
target                  kind             |score|    before     threshold  ratio     ok 
----------------------  ---------------  ---------  ---------  ---------  --------  ---
msm                     fluctuation      0.000e+00  5.344e-03  6.760e-07  0.00e+00  yes
msm[(intercept)]        influence curve  0.000e+00  -          6.760e-07  0.00e+00  yes
msm[assigned contacts]  influence curve  0.000e+00  -          1.569e-07  0.00e+00  yes

PASS: the targeting step solved the estimated efficient score equation.

 bound    psi  delta_from_fitted  ci_lower  ci_upper  truncated_fraction
0.0010 0.2690             0.0002    0.2521    0.2858              0.0003
0.0050 0.2689             0.0001    0.2520    0.2857              0.0023
0.0100 0.2688             0.0000    0.2519    0.2856              0.0100
0.0114 0.2688             0.0000    0.2519    0.2856              0.0120
0.0250 0.2687            -0.0001    0.2519    0.2855              0.0403
0.0500 0.2700             0.0013    0.2536    0.2865              0.1080
0.1000 0.2699             0.0012    0.2546    0.2852              0.2677
0.2000 0.2654            -0.0033    0.2525    0.2784              0.6093
largest absolute slope movement: 0.0033

What this output tells you. Read the parts in order.

output part	what it shows on this draw
the status table and needs attention	one warning, support. The score-equation check passed. A completed row ran, and is not a pass. The deferred row, refute, runs on request
positivity and overlap	36 units (1.20%) sit at the truncation bound. The support report warns when that share is above 1%. The narrowest cadence, high, has an ESS / n of 0.477
the largest clever covariate, 268.5	the high contact count, 6, divided by the smallest fitted probability of high among high patients. It shows that the weight concentrates on a few rows
score-equation check	one row for the fluctuation’s two-coefficient score, and one influence-curve row per coefficient
truncation curve	the slope stays between 0.2654 and 0.2700 across bounds from 0.0010 to 0.2000

The largest bound clips a fraction of 0.6093 of the rows, and the slope moves by at most 0.0033. On this draw, the slope depends little on the truncation bound. That does not verify positivity.

The support report describes the fitted treatment model. It cannot show how small the true probabilities are. About one patient in a hundred has a true probability of high or medium below 0.01. The probe script truth.py gives the true probabilities.

Step 13: sensitivity on the cadence contrasts

Sensitivity analysis asks how strong an unmeasured confounder would need to be to change the conclusion. Sensitivity analysis defines it. cleverly implements the omitted-variable bound for arm-indexed estimands, not for an MSM coefficient. The code reads the robustness value of each contrast of the ATE(reference="low") fit from Step 11. Those contrasts have the same influence curves as the saturated coefficients.

robustness = {
    name: arm_contrasts.sensitivity.robustness_value(estimand=name)
    for name in ("ate[medium vs low]", "ate[high vs low]")
}
for name, values in robustness.items():
    print(
        f"{name}: robustness value {values['rv']:.3f} (confidence-limit value {values['rva']:.3f})"
    )

ate[medium vs low]: robustness value 0.204 (confidence-limit value 0.181)
ate[high vs low]: robustness value 0.420 (confidence-limit value 0.394)

What this output tells you. Each robustness value assumes worst-case alignment (rho=1). It is the equal outcome-side and treatment-side strength that moves the point estimate to zero. The confidence-limit value is the strength at which the confidence limit reaches zero. Both strengths are partial R² values between 0 and 1.

contrast	robustness value	confidence-limit value
ate[medium vs low]	0.204	0.181
ate[high vs low]	0.420	0.394

The medium versus low contrast reaches zero at the weaker confounding strength. These values bound the cadence contrasts, not the slope. The slope is a fixed contrast of the three means (Step 7), but this bound does not cover that combination. The omitted-variable bounds section defines each quantity.

The bound needs a finite second moment of the Riesz representer, so \(E[1/g]\) must be finite for each cadence. On this law it is finite for all three. The review probe script truth.py computes it in closed form from the cadence logits.

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, each cadence’s mean	Ordinary multi-arm TMLE	another law
Step 7, the trend	Point-treatment MSM projection	another law

No registered study covers three cadences, fitted nuisances, or simultaneous bands for a projection.

A review probe is a script under reviews/notebook-review/probes/. The step readings quote its results. A review probe is not a registered study.

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

What the analysis tells the board

question	answer from this analysis
what the slope means	the slope of the least-squares line fitted to the three counterfactual means, with equal weight per cadence: 0.2688 score units per assigned contact
the population and the intervention	adults with a discharge-home order, each assigned one cadence and analyzed as assigned
the uncertainty	the pointwise 95% interval (0.2519, 0.2856), from the influence curve. It is experimental. The review probe measured its coverage at 0.9405 on this law
whether one line is adequate	no. The miss at medium has the interval (-0.2566, -0.1602), which excludes zero
what supports the reading	the score-equation check, the saturated control, and the stable truncation curve
what the diagnostics flagged	the support warning: the fit truncated 1.20% of units, and about 1 patient in 100 has a true probability below 0.01
what supports the cadence contrasts	their robustness values, 0.204 and 0.420. These values do not cover the slope
what remains unverified	no unmeasured confounding, consistency, no interference, positivity beyond the fitted model, the handling of death, and, in a real program, the specification of Q and g

Report the three cadence means or the two contrasts beside the slope. The gain per contact is not constant, and the slope hides that. The analysis does not support these conclusions:

conclusion	why the analysis does not support it
each extra contact raises the score by 0.2688	the line does not fit, so the slope is a summary, not a per-contact effect
ten contacts would raise the score further	no patient was assigned more than six contacts
the board should choose the high cadence	that choice needs the cadence contrasts, their costs, and the board’s own trade-offs
the slope is robust to unmeasured confounding	the omitted-variable bound covers the cadence contrasts, not the slope
the interval covers at 95% in general	one review probe on one law measured 0.9405, and no registered study covers this configuration

Where to go next

The same projection works over regimens and horizons in a longitudinal fit. Its design callable also receives the horizon, and it codes each regimen as a number in the same way as Step 8. A longitudinal MSM projection cross-fits by default, with one pooled update per node. Longitudinal MSM projections shows the call. Sequential offers and time-to-event outcomes estimate the per-regimen parameters that the projection summarizes.

The examples index lists every tutorial in the program.