{ "cells": [ { "cell_type": "markdown", "id": "title", "metadata": {}, "source": [ "# MSM projections: summarize three navigation cadences with one trend\n", "\n", "A marginal structural model (MSM) projection fits a declared working model to counterfactual\n", "means. This notebook fits a straight line in assigned contacts to the means of three navigation\n", "cadences. Each step shows its code, its output, and what the output tells you.\n", "[MSM projections](../technical-reference/msm-projections.md) gives the projection, its clever\n", "covariate, and the algorithm.\n", "\n", "## The situation\n", "\n", "A hospital network sends adults home after a discharge-home order. At that order, staff assign one\n", "of three navigation cadences. Nobody randomized the assignment.\n", "\n", "| tier | assigned protocol | assigned contacts in 30 days |\n", "| --- | --- | --- |\n", "| `low` | one transition-planning contact | 1 |\n", "| `medium` | planning plus one follow-up contact | 2 |\n", "| `high` | planning plus five follow-up contacts | 6 |\n", "\n", "The table separates what the data record from what the question asks about.\n", "\n", "| quantity | kind | when it is fixed | column |\n", "| --- | --- | --- | --- |\n", "| discharge risk, age, baseline support need | observed baseline covariates, W | before assignment | `discharge_risk`, `age`, `baseline_support_need` |\n", "| assigned cadence | observed treatment, A | at the discharge-home order, which is time zero | `cadence` |\n", "| transition score | observed outcome, Y | 30 days after discharge | `transition_score` |\n", "| 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 |\n", "\n", "Completed contacts, readmission, and death occur after assignment. The data set does not record\n", "them, and the analysis must not adjust for them.\n", "\n", "### The causal graph\n", "\n", "```text\n", "discharge_risk ──┬───────────────────► transition_score ◄── baseline_support_need\n", " │ ▲ ▲\n", " ▼ │ │\n", " cadence ───────────────────┘ │\n", " ▲ │\n", "age ─────────────┴────────────────────────────────┘\n", "```\n", "\n", "| arrow | what it asserts | what the analysis does with it |\n", "| --- | --- | --- |\n", "| `discharge_risk` → `cadence`, `age` → `cadence` | staff use risk and age to choose the cadence | with the next row, it makes both variables confounders |\n", "| `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 |\n", "| `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 |\n", "| `cadence` → `transition_score` | the effect that the projection summarizes | the only open path from cadence to score after adjustment |\n", "\n", "| absent arrow | what it asserts | what the analysis does with it |\n", "| --- | --- | --- |\n", "| 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 |\n", "| 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 |\n", "| 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 |\n", "\n", "The graph has no censoring node, because every patient has a score. In this synthetic law, the\n", "cadence effect does not vary with the covariates, so no variable is an effect modifier. A real\n", "program would need to check that claim. The graph says nothing about the shape of the cadence\n", "effect. A straight line in contacts is therefore a choice of summary, not a property of the graph.\n", "\n", "## What you will learn\n", "\n", "| after this notebook you can | the step that shows it |\n", "| --- | --- |\n", "| say why unadjusted cadence means mislead here | Step 3 |\n", "| estimate one counterfactual mean per cadence | Step 6 |\n", "| write the board's question as a projection and fit it | Step 7 |\n", "| code the cadences for the question the board asks | Step 8 |\n", "| test whether the line fits the three means | the failure mode, Step 9 |\n", "| say why the projection weight is part of the estimand | Step 10 |\n", "| check the projection against the per-cadence report | Step 11 |\n", "| read the support, truncation, and sensitivity evidence | Steps 12 and 13 |" ] }, { "cell_type": "markdown", "id": "plan", "metadata": {}, "source": [ "## The causal question and the estimand\n", "\n", "The program board wants one number for its report: how the mean 30-day score changes with assigned\n", "contacts across the three cadences.\n", "\n", "| element | definition |\n", "| --- | --- |\n", "| target population | adults with a discharge-home order at a participating hospital. Step 4 records the eligibility |\n", "| interventions | assign every eligible patient `low`, `medium`, or `high`, analyzed as assigned |\n", "| counterfactual means | E[Y^a] for each cadence a |\n", "| 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 |\n", "\n", "The slope β1 has a value whether or not the three means lie on a line. When they do not, β1 is a\n", "summary of the means. It is not the effect of one more contact.\n", "\n", "| simpler estimand | why it does not answer the board's question |\n", "| --- | --- |\n", "| the difference in observed mean scores | confounding biases it. Step 3 shows the observed order of two cadences reversed |\n", "| 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 |\n", "| 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] |\n", "| 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) |\n", "\n", "Identification needs four assumptions.\n", "\n", "| assumption | what it means for this program | basis |\n", "| --- | --- | --- |\n", "| consistency | each cadence always means its declared script and contacts | scientific judgment |\n", "| no interference | one patient's cadence does not change another patient's score | scientific judgment |\n", "| no unmeasured confounding | risk, age, and support need block every common cause of cadence and score | scientific judgment |\n", "| positivity | each covariate profile has some chance of each cadence | checkable in part, through the support report of Step 12 |\n", "\n", "\n", "\n", "### Comparisons on this page\n", "\n", "Each comparison isolates one reason why two numbers differ. The page names that reason before the\n", "output.\n", "\n", "| step | comparison | estimands | mechanism that makes them differ | evidence |\n", "| --- | --- | --- | --- | --- |\n", "| 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 |\n", "| 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 |\n", "| 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 |\n", "| 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 |\n", "| 11 | the saturated projection against the cadence contrasts | the same estimand, written two ways | none. The two must agree | equal estimates and influence curves |\n", "\n", "## Why this method\n", "\n", "Each method choice answers one feature of the problem.\n", "\n", "| feature of the problem | method choice | reason |\n", "| --- | --- | --- |\n", "| one assignment at time zero, with baseline confounders only | point-treatment TMLE | no later treatment decision, time-varying confounder, or censoring needs longitudinal TMLE |\n", "| a treatment with three unordered labels | a multinomial treatment model g | each cadence gets its own probability, and positivity applies to each cadence |\n", "| 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 |\n", "| 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 |\n", "| a score observed for every patient | no missing-outcome model | no score is missing in this data set |\n", "| finite-dimensional learners and no documented score range | nuisances fitted in sample, without cross-fitting | Step 6 gives the reason |\n", "\n", "The projection has three costs. Its coefficients mean what the working model says. You code each\n", "cadence label as a number yourself. The slope is not the effect of one more contact unless the\n", "three means lie on a line.\n", "\n", "| term | plain meaning |\n", "| --- | --- |\n", "| estimand | the number the question asks for, written before any model is chosen. See [estimands](../user-guide/estimands.md#marginal-structural-models) |\n", "| counterfactual mean | the mean score if every eligible patient were assigned one cadence |\n", "| 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 |\n", "| projection weight | how much each cadence counts when the line is fitted to the means. Here every cadence counts equally |\n", "| 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](../technical-reference/point-treatment-tmle.md) |\n", "| saturated model | a working model with one free coefficient per cadence, so it reproduces the three means exactly |" ] }, { "cell_type": "markdown", "id": "setup-heading", "metadata": {}, "source": [ "## Step 1: set up\n", "\n", "The setup imports the learners and prints the installed `cleverly` version. Every fit below passes\n", "linear learners, a fold count, and a random seed. A rerun therefore reproduces the stored outputs." ] }, { "cell_type": "code", "execution_count": 1, "id": "setup", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:22.055447Z", "iopub.status.busy": "2026-10-16T01:17:22.055137Z", "iopub.status.idle": "2026-10-16T01:17:23.202131Z", "shell.execute_reply": "2026-10-16T01:17:23.201578Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "cleverly 0.1.3\n" ] } ], "source": [ "import numpy as np\n", "import pandas as pd\n", "from sklearn.linear_model import LinearRegression, LogisticRegression\n", "\n", "import cleverly\n", "\n", "pd.set_option(\"display.width\", 130)\n", "pd.set_option(\"display.max_columns\", 20)\n", "print(\"cleverly\", cleverly.__version__)" ] }, { "cell_type": "markdown", "id": "setup-reading", "metadata": {}, "source": [ "**What this output tells you.** The stored outputs in this notebook came from the version named\n", "above. A different version can print different numbers." ] }, { "cell_type": "markdown", "id": "data-heading", "metadata": {}, "source": [ "## Step 2: the data\n", "\n", "`make_multi_arm` draws a three-arm confounded law with known counterfactual means. The code renames\n", "its columns to the program's names. It records the contact count of each cadence, prints the first\n", "rows, and prints the known means." ] }, { "cell_type": "code", "execution_count": 2, "id": "data", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.203898Z", "iopub.status.busy": "2026-10-16T01:17:23.203260Z", "iopub.status.idle": "2026-10-16T01:17:23.436067Z", "shell.execute_reply": "2026-10-16T01:17:23.435202Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "rows and columns: (3000, 5)\n", "cadence labels: ['high', 'low', 'medium']\n", " transition_score cadence discharge_risk age baseline_support_need\n", "0 1.266 high -0.567 -1.057 1.062\n", "1 -0.691 low -1.129 -1.779 0.655\n", "2 1.325 medium 0.841 -0.949 -0.565\n", "3 5.003 medium 2.780 -1.017 -0.506\n", "4 0.847 medium 1.179 0.745 1.612\n", "\n", "known counterfactual means of the synthetic law:\n", " ey[low] (1 contacts): -0.0000\n", " ey[medium] (2 contacts): 0.6000\n", " ey[high] (6 contacts): 1.4400\n", "gain per contact, low to medium: 0.60\n", "gain per contact, medium to high: 0.21\n" ] } ], "source": [ "from cleverly.datasets import make_multi_arm\n", "\n", "frame, truth = make_multi_arm(n=3_000, seed=61)\n", "frame = frame.rename(\n", " columns={\n", " \"Y\": \"transition_score\",\n", " \"A\": \"cadence\",\n", " \"W1\": \"discharge_risk\",\n", " \"W2\": \"age\",\n", " \"W3\": \"baseline_support_need\",\n", " }\n", ")\n", "ARMS = (\"low\", \"medium\", \"high\")\n", "CONTACTS_30D = {\"low\": 1.0, \"medium\": 2.0, \"high\": 6.0}\n", "contacts = np.array([CONTACTS_30D[arm] for arm in ARMS])\n", "population = np.array([truth[f\"ey[{arm}]\"] for arm in ARMS])\n", "\n", "print(\"rows and columns:\", frame.shape)\n", "print(\"cadence labels:\", sorted(frame[\"cadence\"].unique()))\n", "print(frame.head().round(3))\n", "print()\n", "print(\"known counterfactual means of the synthetic law:\")\n", "for arm, value in zip(ARMS, population, strict=True):\n", " print(f\" ey[{arm}] ({CONTACTS_30D[arm]:.0f} contacts): {value:.4f}\")\n", "gain_per_contact = np.diff(population) / np.diff(contacts)\n", "print(f\"gain per contact, low to medium: {gain_per_contact[0]:.2f}\")\n", "print(f\"gain per contact, medium to high: {gain_per_contact[1]:.2f}\")" ] }, { "cell_type": "markdown", "id": "data-reading", "metadata": {}, "source": [ "**What this output tells you.** Each row is one discharge, and `cadence` holds a text label. The\n", "published truths use the generator's labels, so the labels stay unchanged.\n", "\n", "The generator draws the three covariates as independent standard-normal variables. Risk and age\n", "enter both the cadence probabilities and the score. Support need enters only the score. These are\n", "the arrows of the causal graph. The score has no unit and no documented range.\n", "\n", "The true means are 0 for `low` (printed as `-0.0000` after rounding), 0.6000 for `medium`, and\n", "1.4400 for `high`. The gain per contact falls from 0.60 between `low` and `medium` to 0.21 between\n", "`medium` and `high`. No straight line in assigned contacts passes through the three means, so the\n", "working model of Step 7 is misspecified on purpose.\n", "\n", "A real program has no `truth`. Every comparison against it below is a teaching device." ] }, { "cell_type": "markdown", "id": "association-heading", "metadata": {}, "source": [ "## Step 3: association first\n", "\n", "A confounder is a variable that changes both the assigned cadence and the outcome. The code\n", "compares the cadences before any adjustment. It prints the observed mean score, the means of the\n", "two confounders, the share of patients, and the true counterfactual mean." ] }, { "cell_type": "code", "execution_count": 3, "id": "association", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.437758Z", "iopub.status.busy": "2026-10-16T01:17:23.437496Z", "iopub.status.idle": "2026-10-16T01:17:23.451114Z", "shell.execute_reply": "2026-10-16T01:17:23.450391Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " transition_score discharge_risk age share true mean\n", "cadence \n", "low -0.082 -0.132 -0.128 0.272 -0.00\n", "medium 1.324 0.532 -0.343 0.368 0.60\n", "high 0.808 -0.387 0.474 0.360 1.44\n" ] } ], "source": [ "observed = frame.groupby(\"cadence\")[[\"transition_score\", \"discharge_risk\", \"age\"]].mean()\n", "observed = observed.loc[list(ARMS)]\n", "observed[\"share\"] = frame[\"cadence\"].value_counts(normalize=True).loc[list(ARMS)]\n", "observed[\"true mean\"] = population\n", "print(observed.round(3))" ] }, { "cell_type": "markdown", "id": "association-reading", "metadata": {}, "source": [ "**What this output tells you.** The observed means put `medium` above `high`, at 1.324 against\n", "0.808. The true means reverse that order, at 0.60 against 1.44.\n", "\n", "The confounded assignment causes the reversal. The table compares the two confounders in the `medium` and\n", "`high` groups.\n", "\n", "| confounder | `medium` mean | `high` mean | effect on the score in this law |\n", "| --- | --- | --- | --- |\n", "| `discharge_risk` | 0.532 | -0.387 | a higher value raises the score |\n", "| `age` | -0.343 | 0.474 | a higher value lowers the score |\n", "\n", "Both confounders push the observed `medium` mean above the observed `high` mean. A trend fitted to\n", "these observed means would summarize the confounding too.\n", "\n", "The shares are 0.272, 0.368, and 0.360. Step 10 uses them to show that the projection weight is part\n", "of the estimand." ] }, { "cell_type": "markdown", "id": "protocol-heading", "metadata": {}, "source": [ "## Step 4: write the protocol\n", "\n", "A `StudyProtocol` records the scientific design before any model runs. The record gets a\n", "fingerprint, and every result fitted from it carries that fingerprint. This page starts from\n", "`navigation_protocol()`, the protocol of the\n", "[shared study design](index.md#the-shared-study-design). `dataclasses.replace` lists the three\n", "cadences as three strategies, each with its version." ] }, { "cell_type": "code", "execution_count": 4, "id": "protocol", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.453027Z", "iopub.status.busy": "2026-10-16T01:17:23.452585Z", "iopub.status.idle": "2026-10-16T01:17:23.457779Z", "shell.execute_reply": "2026-10-16T01:17:23.457054Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "causal study protocol: schema 1; ec77f97f713c2306\n", "target population: Adults with a discharge-home order at a participating hospital during the enrollment period\n", "eligibility: ['Age 18 years or older', 'Discharge home ordered at a participating hospital']\n", "time zero: Discharge-home order, after baseline measurement and before cadence assignment\n", "treatment strategies: ['Assign the low navigation cadence', 'Assign the medium navigation cadence', 'Assign the high navigation cadence']\n", "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']\n", "outcome: Patient-reported transition score in synthetic units, with no documented range\n", "horizon: 30 days after discharge\n", "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']\n", "interference unit: Individual patient\n", "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']\n" ] } ], "source": [ "from dataclasses import replace\n", "\n", "from cleverly.datasets import navigation_protocol\n", "\n", "program = navigation_protocol()\n", "protocol = replace(\n", " program,\n", " outcome=\"Patient-reported transition score in synthetic units, with no documented range\",\n", " time_zero=\"Discharge-home order, after baseline measurement and before cadence assignment\",\n", " treatment_strategies=(\n", " \"Assign the low navigation cadence\",\n", " \"Assign the medium navigation cadence\",\n", " \"Assign the high navigation cadence\",\n", " ),\n", " treatment_versions=(\n", " \"One transition-planning contact within 30 days\",\n", " \"Transition planning plus one follow-up contact within 30 days\",\n", " \"Transition planning plus five follow-up contacts within 30 days\",\n", " ),\n", " intercurrent_event_handling=(\n", " program.intercurrent_event_handling[0],\n", " \"Analyze the assigned cadence regardless of completed contacts\",\n", " \"The score has no worst value, so the protocol handles death before day 30 \"\n", " \"by the hypothetical strategy: the score had the patient not died\",\n", " ),\n", " assumption_rationale=(\n", " \"The recorded baseline variables cover the measured common causes of cadence and score\",\n", " \"Every cadence uses the same script and access rules, which supports consistency\",\n", " program.assumption_rationale[2],\n", " ),\n", ")\n", "print(\"\\n\".join(protocol.summary_lines()))" ] }, { "cell_type": "markdown", "id": "protocol-reading", "metadata": {}, "source": [ "**What this output tells you.** The first line gives the schema version and the fingerprint\n", "`ec77f97f713c2306`. The other lines repeat each field. This page changes six fields of\n", "`navigation_protocol()`, and\n", "[point-treatment TMLE](point-treatment-tmle.ipynb#step-4-write-the-protocol) reads the others.\n", "\n", "| protocol field | what this page writes |\n", "| --- | --- |\n", "| outcome | the score is in synthetic units and has no documented range. The program records a share of the maximum score |\n", "| time zero | the discharge-home order, before cadence assignment |\n", "| treatment strategies and versions | the three cadences, with their contacts in words |\n", "| 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 |\n", "| assumption rationale | the common causes of cadence and score, and one script for every cadence |\n", "\n", "Under the hypothetical strategy, a real analysis treats death as censoring. It then needs\n", "exchangeability and positivity for death given the recorded history. This synthetic law draws no\n", "deaths, so the choice does not change the estimate here.\n", "\n", "A score with no documented range has no range to declare as `q_bounds`. Step 6 uses that fact.\n", "\n", "The protocol has no field for the contact mapping, the working model, or the projection weight. The\n", "typed estimand in Step 7 owns all three." ] }, { "cell_type": "markdown", "id": "identify-heading", "metadata": {}, "source": [ "## Step 5: design and identification\n", "\n", "The design is an ordinary point-treatment design with a three-level treatment. `CounterfactualMean`\n", "asks for one mean per cadence. The projection in Step 7 summarizes those means." ] }, { "cell_type": "code", "execution_count": 5, "id": "identify", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.459686Z", "iopub.status.busy": "2026-10-16T01:17:23.459485Z", "iopub.status.idle": "2026-10-16T01:17:23.474336Z", "shell.execute_reply": "2026-10-16T01:17:23.473599Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "counterfactual mean under each treatment, E[Y^a]\n", "identified by explicit-adjustment: E_W[E(transition_score | cadence=a, W)] for a in ['high', 'low', 'medium']\n", "adjustment/history: ['discharge_risk', 'age', 'baseline_support_need']\n", "required nuisances: ['outcome_regression', 'treatment_mechanism']\n", "assumptions:\n", " - consistency: Y = Y^a when A = a\n", " - no interference: one unit's potential outcome does not depend on other units' treatment assignments\n", " - no unmeasured confounding: Y^a is independent of A given W\n", " - positivity: P(cadence = a | W) > 0 almost surely for every supported treatment level a in ['high', 'low', 'medium']\n", "causal study protocol: schema 1; ec77f97f713c2306\n" ] } ], "source": [ "from cleverly import CausalStudy, CounterfactualMean, PointTreatment\n", "\n", "study = CausalStudy(\n", " frame,\n", " design=PointTreatment(\n", " outcome=\"transition_score\",\n", " treatment=\"cadence\",\n", " adjustment=(\"discharge_risk\", \"age\", \"baseline_support_need\"),\n", " ),\n", " protocol=protocol,\n", ")\n", "arm_means_effect = study.identify(CounterfactualMean())\n", "print(arm_means_effect.summary(protocol=\"fingerprint\"))" ] }, { "cell_type": "markdown", "id": "identify-reading", "metadata": {}, "source": [ "**What this output tells you.** The formula averages the outcome regression over all patients, once\n", "for each cadence. The required nuisances are the outcome regression and the treatment mechanism.\n", "The four assumptions are the ones that\n", "[the causal question](#the-causal-question-and-the-estimand) lists with their basis.\n", "\n", "The positivity condition names `high`, `low`, and `medium`. Each level needs a positive probability\n", "at every covariate value where the target places mass. The summary then prints the fingerprint of\n", "the stored protocol." ] }, { "cell_type": "markdown", "id": "estimate-arms-heading", "metadata": {}, "source": [ "## Step 6: estimate each cadence's mean\n", "\n", "Start with the per-cadence report, because the projection is a summary of it. The linear outcome\n", "model and the multinomial logistic treatment model are correctly specified for this law.\n", "\n", "| setting | value | what it does |\n", "| --- | --- | --- |\n", "| `outcome_learner` | `LinearRegression()` | fits Q, the expected score given the cadence and the covariates |\n", "| `treatment_learner` | `LogisticRegression(max_iter=1000, random_state=61)` | fits g, the probability of each cadence given the covariates |\n", "| `CrossFitting(enabled=False)` | no outer folds | fits each nuisance on every row |\n", "| `Runtime(random_state=61, n_jobs=1)` | fixed seed, one process | makes the fit reproducible |\n", "\n", "The fit runs in sample for two reasons. A cross-fitted fit of a continuous outcome declares the\n", "outcome's range as `q_bounds`, and this score has no documented range. The learners are\n", "finite-dimensional regressions. That makes plausible the data-reuse condition that\n", "[CV-TMLE](../technical-reference/cv-tmle.md#what-this-solves) removes. The\n", "[cross-fitting tutorial](cross-fitting.ipynb) shows the fold layer on a score with a recorded range.\n", "\n", "Targeting makes a small update to Q, weighted by g, that removes first-order bias.\n", "[Targeting and bounds](../user-guide/methods-learners.md#targeting-and-bounds) defines it." ] }, { "cell_type": "code", "execution_count": 6, "id": "estimate-arms", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.476955Z", "iopub.status.busy": "2026-10-16T01:17:23.476560Z", "iopub.status.idle": "2026-10-16T01:17:23.583194Z", "shell.execute_reply": "2026-10-16T01:17:23.582905Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " psi std_err ci_lower ci_upper true mean\n", "estimand \n", "ey[low] -0.0038 0.0365 -0.0754 0.0678 -0.00\n", "ey[medium] 0.6399 0.0386 0.5641 0.7156 0.60\n", "ey[high] 1.4642 0.0390 1.3877 1.5408 1.44\n", "\n", "protocol fingerprint on the fit: ec77f97f713c2306\n" ] } ], "source": [ "from cleverly import CrossFitting, ModelSpec, Runtime, TMLEMethod\n", "\n", "method = TMLEMethod(\n", " models=ModelSpec(\n", " outcome_learner=LinearRegression(),\n", " treatment_learner=LogisticRegression(max_iter=1000, random_state=61),\n", " ),\n", " cross_fitting=CrossFitting(enabled=False),\n", " runtime=Runtime(random_state=61, n_jobs=1),\n", ")\n", "arm_result = arm_means_effect.estimate(method=method)\n", "\n", "arm_table = arm_result.to_frame().set_index(\"estimand\").loc[[f\"ey[{arm}]\" for arm in ARMS]]\n", "arm_table = arm_table[[\"psi\", \"std_err\", \"ci_lower\", \"ci_upper\"]]\n", "arm_table[\"true mean\"] = population\n", "print(arm_table.round(4))\n", "print()\n", "print(\"protocol fingerprint on the fit:\", arm_result.provenance.protocol_fingerprint)" ] }, { "cell_type": "markdown", "id": "estimate-arms-reading", "metadata": {}, "source": [ "**What this output tells you.** Each row is one counterfactual mean. The `medium` estimate is 0.6399\n", "with the interval (0.5641, 0.7156), and the true mean is 0.60. The `high` estimate is 1.4642, and\n", "the `low` estimate is -0.0038.\n", "\n", "Each interval contains its true mean on this draw. That is one draw, not a coverage result. The\n", "interval comes from the influence curve, which measures how much each row moves the estimate. See\n", "[inference](../technical-reference/inference.md). The fit carries the protocol fingerprint\n", "`ec77f97f713c2306`.\n" ] }, { "cell_type": "markdown", "id": "trend-heading", "metadata": {}, "source": [ "## Step 7: declare the working model and fit the trend\n", "\n", "The `design` function receives one cadence label and the covariate frame. It returns the design\n", "matrix for that cadence, with one column per term. The `terms` tuple names the columns.\n", "\n", "The declaration leaves `weights` unset, which sets the projection weight h(a, V) to 1 for every\n", "cadence. The same `method` fits the projection." ] }, { "cell_type": "code", "execution_count": 7, "id": "trend", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.596200Z", "iopub.status.busy": "2026-10-16T01:17:23.596013Z", "iopub.status.idle": "2026-10-16T01:17:23.670337Z", "shell.execute_reply": "2026-10-16T01:17:23.670029Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Targeted maximum likelihood estimation\n", "======================================\n", "n = 3000; covariates = 3; arm shares: high=0.36, low=0.272, medium=0.368\n", "causal estimand: projection of counterfactual means onto a working model\n", "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\n", "required nuisances: outcome_regression, treatment_mechanism\n", "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\n", "causal study protocol: schema 1; ec77f97f713c2306\n", "outcome family: gaussian; fluctuation: logistic (iterative)\n", "TMLE (in-sample nuisances): nuisances fitted in-sample (cross_fit=False)\n", "propensity truncated to [0.0114, 0.9886]; ATT/ATC to [0.025, 0.975]\n", "outcome scaled from [-6.044, 7.723] to [0, 1]\n", "data 248bf6e493a8a2be | folds 9ba86d4d13cb065e | cleverly 0.1.3 | 2026-10-16T01:17:23+00:00\n", "\n", "estimand psi std_err 95% CI p_value support \n", "---------------------- ------- -------- -------------------- ------- ------------\n", "msm[(intercept)] -0.1068 0.03703 [-0.17937, -0.03423] 0.0039 experimental\n", "msm[assigned contacts] 0.26875 0.008594 [0.25191, 0.2856] <1e-4 experimental\n", "\n", "experimental: No registered study covers this configuration. Nearest row: T18 (point-msm), differing on treatment.\n", "\n", "simultaneous 95% bands (multiplier bootstrap, critical value 2.107 vs 1.960 pointwise) [experimental]:\n", " msm[(intercept)] [-0.18484, -0.028769]\n", " msm[assigned contacts] [0.25064, 0.28686]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "from cleverly import MSMProjection\n", "from cleverly.msm import MSM\n", "\n", "\n", "def contacts_design(arm, data):\n", " \"\"\"One row per patient: an intercept and the cadence's assigned contact count.\"\"\"\n", " return np.column_stack([np.ones(len(data)), np.full(len(data), CONTACTS_30D[arm])])\n", "\n", "\n", "trend = MSM(design=contacts_design, terms=(\"(intercept)\", \"assigned contacts\"), design_kind=\"known\")\n", "trend_effect = study.identify(MSMProjection(trend))\n", "trend_result = trend_effect.estimate(method=method)\n", "print(trend_result.summary(protocol=\"fingerprint\"))" ] }, { "cell_type": "markdown", "id": "trend-reading", "metadata": {}, "source": [ "**What this output tells you.** The `identification assumptions` line lists the four causal\n", "assumptions from Step 5 and adds two.\n", "\n", "| added assumption | what it means here |\n", "| --- | --- |\n", "| 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 |\n", "| the weighted Gram matrix is invertible | the three cadences have distinct contact counts, so one line is closest |\n", "\n", "The configuration lines have the same meaning as in\n", "[point-treatment TMLE](point-treatment-tmle.ipynb#step-6-estimate-the-ate). Here the fluctuation\n", "solves one score equation per coefficient, not one per cadence.\n", "\n", "The slope is 0.26875 transition-score units per additional assigned contact. The intercept is\n", "-0.1068, the line's value at zero contacts, which no cadence assigns.\n", "\n", "The simultaneous 95% bands are built to contain both coefficients at the same time, so their\n", "critical value is 2.107 instead of 1.960. Quote the pointwise interval when you report the slope\n", "alone.\n", "\n", "The `support` column marks both coefficients and the bands `experimental`, because no registered\n", "study covers this configuration. The [support label](../user-guide/results-assessment.md#read-the-support-label) makes no coverage claim. The review probe\n", "`reviews/notebook-review/probes/msm-projections-final/sweep.py` repeats this fit on 4,000 draws of\n", "this law at n = 3,000. The interval covered on 3,762 of 4,000 draws, a rate of 0.9405. The probe\n", "does not isolate the cause of the shortfall.\n", "\n", "With a uniform weight, the slope is a fixed contrast of the three means:\n", "(-2 E[Y(low)] - E[Y(medium)] + 3 E[Y(high)]) / 14. Step 9 checks that identity on the population\n", "means." ] }, { "cell_type": "markdown", "id": "msm-step-coding-heading", "metadata": {}, "source": [ "## Step 8: code the cadences for the question the board asks\n", "\n", "The contact counts are the program's decision, not the estimator's. The number that a `design`\n", "function writes for each label sets the question the slope answers. The code fits a second known\n", "design that codes each cadence by its step, 0, 1, and 2. It prints both slopes beside the population\n", "slope of each coding." ] }, { "cell_type": "code", "execution_count": 8, "id": "msm-step-coding", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.677931Z", "iopub.status.busy": "2026-10-16T01:17:23.677777Z", "iopub.status.idle": "2026-10-16T01:17:23.747337Z", "shell.execute_reply": "2026-10-16T01:17:23.746784Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coefficient low, medium, high slope std_err 95% CI population slope\n", "msm[assigned contacts] 1, 2, 6 0.2688 0.0086 (0.2519, 0.2856) 0.2657\n", " msm[cadence step] 0, 1, 2 0.7341 0.0224 (0.6902, 0.7779) 0.7200\n" ] } ], "source": [ "CADENCE_STEP = {\"low\": 0.0, \"medium\": 1.0, \"high\": 2.0}\n", "\n", "\n", "def cadence_step_design(arm, data):\n", " \"\"\"One row per patient: an intercept and the cadence's step, 0, 1, or 2.\"\"\"\n", " return np.column_stack([np.ones(len(data)), np.full(len(data), CADENCE_STEP[arm])])\n", "\n", "\n", "step_trend = MSM(\n", " design=cadence_step_design, terms=(\"(intercept)\", \"cadence step\"), design_kind=\"known\"\n", ")\n", "step_result = study.identify(MSMProjection(step_trend)).estimate(method=method)\n", "\n", "codings = (\n", " (\"assigned contacts\", CONTACTS_30D, trend_result),\n", " (\"cadence step\", CADENCE_STEP, step_result),\n", ")\n", "coding_rows = []\n", "for term, coding, result in codings:\n", " values = np.array([coding[arm] for arm in ARMS])\n", " coding_line, *_ = np.linalg.lstsq(np.column_stack([np.ones(3), values]), population, rcond=None)\n", " slope = result[f\"msm[{term}]\"]\n", " coding_rows.append(\n", " {\n", " \"coefficient\": f\"msm[{term}]\",\n", " \"low, medium, high\": \", \".join(f\"{value:.0f}\" for value in values),\n", " \"slope\": f\"{slope.psi:.4f}\",\n", " \"std_err\": f\"{slope.std_error:.4f}\",\n", " \"95% CI\": f\"({slope.ci[0]:.4f}, {slope.ci[1]:.4f})\",\n", " \"population slope\": f\"{coding_line[1]:.4f}\",\n", " }\n", " )\n", "print(pd.DataFrame(coding_rows).to_string(index=False))" ] }, { "cell_type": "markdown", "id": "msm-step-coding-reading", "metadata": {}, "source": [ "**What this output tells you.** The two rows read the same three means through two codings.\n", "\n", "| coefficient | slope | interval | population slope | question it answers |\n", "| --- | --- | --- | --- | --- |\n", "| `msm[assigned contacts]` | 0.2688 | (0.2519, 0.2856) | 0.2657 | the change per assigned contact |\n", "| `msm[cadence step]` | 0.7341 | (0.6902, 0.7779) | 0.7200 | the change per cadence step |\n", "\n", "The step coding counts the move from `medium` to `high` as one unit. The contact coding counts the\n", "same move as four contacts. Neither line passes through the three means. The board asks for the\n", "change per contact, so the rest of this page reads the contact slope. Each interval contains its\n", "population slope on this draw.\n", "\n", "Under the uniform weight, the per-step slope is half the difference between the `high` and `low`\n", "means. Its population value is 1.44 / 2 = 0.7200. Step 11 prints the fitted `high vs low` contrast\n", "as 1.4681, and half of that is 0.734.\n", "\n", "Write the coding in the `design` function, as both fits here do. The labels have no order of their\n", "own. Sorted as text, they are `high`, `low`, and `medium`." ] }, { "cell_type": "markdown", "id": "failure-heading", "metadata": {}, "source": [ "## Step 9: the failure mode, a projection is not a fitted curve\n", "\n", "The working model from Step 7 is wrong for this law. This step compares two estimands: the\n", "projection's value at each cadence, β0 + β1 c(a), and the cadence mean E[Y^a]. The mechanism that\n", "separates them is the falling gain per contact that Step 2 prints. The code compares the line with the means,\n", "first in the population and then in the fits. It measures the miss at `medium` with\n", "`arm_result.contrast(...)`, which reads the joint influence curve of the three cadence means. It\n", "also computes the population slope as the fixed contrast of the three means." ] }, { "cell_type": "code", "execution_count": 9, "id": "failure-mode", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.749684Z", "iopub.status.busy": "2026-10-16T01:17:23.749536Z", "iopub.status.idle": "2026-10-16T01:17:23.785640Z", "shell.execute_reply": "2026-10-16T01:17:23.775818Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "in the population:\n", " contacts population mean population line line minus mean\n", "low 1.0 -0.00 0.1486 0.1486\n", "medium 2.0 0.60 0.4143 -0.1857\n", "high 6.0 1.44 1.4771 0.0371\n", "\n", "in the fits:\n", " estimated mean ci_lower ci_upper estimated line line minus mean\n", "low -0.0038 -0.0754 0.0678 0.1620 0.1658\n", "medium 0.6399 0.5641 0.7156 0.4307 -0.2092\n", "high 1.4642 1.3877 1.5408 1.5057 0.0415\n", "absolute Step 7 miss at medium: 0.2092\n", "\n", "medium: line minus mean, -0.2084 CI=(-0.2566, -0.1602)\n", "absolute difference from the Step 7 miss: 0.0008\n", "\n", "population projection: intercept -0.1171, slope 0.2657\n", "slope as the fixed contrast of the three means: 0.2657\n", "msm[(intercept)]: CI=(-0.1794, -0.0342) population projection -0.1171\n", "msm[assigned contacts]: CI=(0.2519, 0.2856) population projection 0.2657\n" ] } ], "source": [ "design = np.vstack([contacts_design(arm, frame.iloc[:1]) for arm in ARMS])\n", "names = [f\"msm[{term}]\" for term in trend.terms]\n", "projection, *_ = np.linalg.lstsq(design, population, rcond=None)\n", "estimated_line = design @ np.array([trend_result[name].psi for name in names])\n", "\n", "in_population = pd.DataFrame(\n", " {\n", " \"contacts\": contacts,\n", " \"population mean\": population,\n", " \"population line\": design @ projection,\n", " \"line minus mean\": design @ projection - population,\n", " },\n", " index=list(ARMS),\n", ")\n", "in_fits = (\n", " arm_table[[\"psi\", \"ci_lower\", \"ci_upper\"]]\n", " .set_axis(list(ARMS))\n", " .rename(columns={\"psi\": \"estimated mean\"})\n", ")\n", "in_fits[\"estimated line\"] = estimated_line\n", "in_fits[\"line minus mean\"] = estimated_line - in_fits[\"estimated mean\"]\n", "print(\"in the population:\")\n", "print(in_population.round(4).to_string())\n", "print()\n", "print(\"in the fits:\")\n", "print(in_fits.round(4).to_string())\n", "step7_medium_miss = in_fits.loc[\"medium\", \"line minus mean\"]\n", "print(f\"absolute Step 7 miss at medium: {abs(step7_medium_miss):.4f}\")\n", "print()\n", "\n", "# The least-squares line at medium is a fixed weighted sum of the three means.\n", "misfit_weights = (design @ np.linalg.pinv(design))[1] - np.eye(3)[1]\n", "medium_misfit = arm_result.contrast(\n", " lambda means: float(misfit_weights @ means),\n", " [f\"ey[{arm}]\" for arm in ARMS],\n", " name=\"medium: line minus mean\",\n", " gradient=lambda means: misfit_weights,\n", ")\n", "low, high = medium_misfit.ci\n", "print(f\"medium: line minus mean, {medium_misfit.psi:.4f} CI=({low:.4f}, {high:.4f})\")\n", "print(f\"absolute difference from the Step 7 miss: {abs(medium_misfit.psi - step7_medium_miss):.4f}\")\n", "print()\n", "print(f\"population projection: intercept {projection[0]:.4f}, slope {projection[1]:.4f}\")\n", "contrast_slope = (-2 * population[0] - population[1] + 3 * population[2]) / 14\n", "print(f\"slope as the fixed contrast of the three means: {contrast_slope:.4f}\")\n", "for name, target in zip(names, projection, strict=True):\n", " low, high = trend_result[name].ci\n", " print(f\"{name}: CI=({low:.4f}, {high:.4f}) population projection {target:.4f}\")" ] }, { "cell_type": "markdown", "id": "failure-reading", "metadata": {}, "source": [ "**What this output tells you.** In the population, the line misses `medium` by about 0.19, since\n", "`line minus mean` is -0.1857. It overshoots `low` by 0.1486 and `high` by 0.0371.\n", "\n", "In the fits, the `absolute Step 7 miss at medium` line reports 0.2092. The\n", "`medium: line minus mean` row measures that miss with an interval.\n", "\n", "| part of the row | value on this draw | what it means |\n", "| --- | --- | --- |\n", "| estimate | -0.2084 | the line fitted to the three estimated means, minus the estimated `medium` mean |\n", "| interval | (-0.2566, -0.1602) | it excludes zero, so the data reject, at the 5% level, a line through all three means |\n", "| population miss | -0.1857 | the interval contains it |\n", "\n", "With three cadences and two coefficients, a nonzero miss at one cadence is equivalent to a\n", "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\n", "means, so their correlation enters the interval. Its estimate differs from the Step 7 miss by 0.0008\n", "on this draw, because the cadence fit and the MSM fit target separately.\n", "\n", "| check | value on this draw |\n", "| --- | --- |\n", "| population projection | intercept -0.1171, slope 0.2657 |\n", "| slope as the fixed contrast of the means | 0.2657 |\n", "| slope interval from the fit | (0.2519, 0.2856), which contains 0.2657 |\n", "| intercept interval from the fit | (-0.1794, -0.0342), which contains -0.1171 |\n", "\n", "Read the slope as \"the best linear summary of the cadence response under a uniform weight\". Do not\n", "read it as \"the causal effect of one more contact\". No patient was assigned more than six contacts,\n", "so the projection says nothing about a cadence with ten." ] }, { "cell_type": "markdown", "id": "share-weight-heading", "metadata": {}, "source": [ "## Step 10: the weight is part of the estimand\n", "\n", "The projection weight sets how much each cadence counts when the line is fitted. The code fits the\n", "population line again, with each cadence weighted by its share from Step 3. It then declares a fixed\n", "weight of 1, 10, and 1 for `low`, `medium`, and `high`. It prints the population slope under that\n", "weight, and it fits the projection with the weight declared as `weights_kind=\"known\"`." ] }, { "cell_type": "code", "execution_count": 10, "id": "share-weight", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.787719Z", "iopub.status.busy": "2026-10-16T01:17:23.787553Z", "iopub.status.idle": "2026-10-16T01:17:23.879374Z", "shell.execute_reply": "2026-10-16T01:17:23.879008Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "cadence shares: low 0.272, medium 0.368, high 0.360\n", "population slope, uniform weight: 0.2657\n", "population slope, share weight: 0.2593\n", "shift from the share weight: 0.0065\n", "standard error of the Step 7 slope: 0.0086\n", "\n", "population slope, 1:10:1 weight: 0.2400\n", "fitted slope, uniform weight: 0.2688 CI=(0.2519, 0.2856)\n", "fitted slope, 1:10:1 weight: 0.2398 CI=(0.2214, 0.2582)\n" ] } ], "source": [ "shares = observed[\"share\"].to_numpy()\n", "# Weighted least squares: solve (X' W X) b = X' W y, with W the diagonal matrix of weights.\n", "share_projection = np.linalg.solve(\n", " design.T @ (shares[:, None] * design), design.T @ (shares * population)\n", ")\n", "slope = trend_result[\"msm[assigned contacts]\"]\n", "print(\n", " \"cadence shares:\",\n", " \", \".join(f\"{arm} {share:.3f}\" for arm, share in zip(ARMS, shares, strict=True)),\n", ")\n", "print(f\"population slope, uniform weight: {projection[1]:.4f}\")\n", "print(f\"population slope, share weight: {share_projection[1]:.4f}\")\n", "print(f\"shift from the share weight: {abs(share_projection[1] - projection[1]):.4f}\")\n", "print(f\"standard error of the Step 7 slope: {slope.std_error:.4f}\")\n", "print()\n", "\n", "# A fixed weight that puts most of the projection on medium.\n", "FIXED_WEIGHT = {\"low\": 1.0, \"medium\": 10.0, \"high\": 1.0}\n", "fixed = np.array([FIXED_WEIGHT[arm] for arm in ARMS])\n", "fixed_projection = np.linalg.solve(\n", " design.T @ (fixed[:, None] * design), design.T @ (fixed * population)\n", ")\n", "weighted_trend = MSM(\n", " design=contacts_design,\n", " terms=trend.terms,\n", " weights=lambda arm, data: np.full(len(data), FIXED_WEIGHT[arm]),\n", " weights_kind=\"known\",\n", " design_kind=\"known\",\n", ")\n", "weighted_result = study.identify(MSMProjection(weighted_trend)).estimate(method=method)\n", "weighted_slope = weighted_result[\"msm[assigned contacts]\"]\n", "print(f\"population slope, 1:10:1 weight: {fixed_projection[1]:.4f}\")\n", "for label, estimate in ((\"uniform\", slope), (\"1:10:1\", weighted_slope)):\n", " low, high = estimate.ci\n", " print(f\"fitted slope, {label} weight: {estimate.psi:.4f} CI=({low:.4f}, {high:.4f})\")" ] }, { "cell_type": "markdown", "id": "share-weight-reading", "metadata": {}, "source": [ "**What this output tells you.** The uniform weight gives a population slope of 0.2657. The share\n", "weight gives 0.2593. The means are the same, so only the weight moved the slope.\n", "\n", "The shift is 0.0065, and the standard error of the Step 7 slope is 0.0086. That standard error\n", "describes one slope, not the difference between two slopes. This illustration uses population means\n", "and sample shares. It does not test whether two fitted slopes differ.\n", "\n", "The 1:10:1 weight moves the population slope to 0.2400. The fitted slope under that weight is\n", "0.2398, with the interval (0.2214, 0.2582). On this draw, each interval contains its own population\n", "slope and not the other.\n", "\n", "| fitted slope | interval | contains 0.2657 | contains 0.2400 |\n", "| --- | --- | --- | --- |\n", "| uniform weight, 0.2688 | (0.2519, 0.2856) | yes | no |\n", "| 1:10:1 weight, 0.2398 | (0.2214, 0.2582) | no | yes |\n", "\n", "A change of weight or terms changes the estimand, not only the estimate. A misspecified working\n", "model still defines a parameter. Interval validity still needs the causal and nuisance conditions.\n", "\n", "Declare a `weights=` function as known only when it does not depend on the data. The 1:10:1 weight\n", "qualifies. The sample shares and a function of the fitted treatment mechanism do not. Such a weight,\n", "declared as known, still fits. That fit omits the weight's influence-curve term, so its standard\n", "error can be too small or too large. The\n", "[technical entry](../technical-reference/msm-projections.md#variations) names the term that an\n", "estimated weight or design adds." ] }, { "cell_type": "markdown", "id": "saturated-heading", "metadata": {}, "source": [ "## Step 11: the control, a saturated working model\n", "\n", "A saturated model has one free coefficient per cadence. It cannot be misspecified, so it represents\n", "the three means exactly. The code fits it beside the average treatment effect (`ATE`) contrasts\n", "against `low`, with the same method. Each row pairs one coefficient with its counterpart." ] }, { "cell_type": "code", "execution_count": 11, "id": "saturated", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:23.888586Z", "iopub.status.busy": "2026-10-16T01:17:23.888304Z", "iopub.status.idle": "2026-10-16T01:17:24.023418Z", "shell.execute_reply": "2026-10-16T01:17:24.022980Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coefficient psi ci counterpart counterpart psi counterpart ci influence curves agree within 1e-12\n", " msm[(intercept)] -0.0038 (-0.0754, 0.0678) ey[low] -0.0038 (-0.0754, 0.0678) True\n", "msm[medium vs low] 0.6437 (0.5562, 0.7312) ate[medium vs low] 0.6437 (0.5562, 0.7312) True\n", " msm[high vs low] 1.4681 (1.3805, 1.5557) ate[high vs low] 1.4681 (1.3805, 1.5557) True\n" ] } ], "source": [ "from cleverly import ATE\n", "\n", "\n", "def cadence_indicators(arm, data):\n", " \"\"\"One free coefficient per cadence: an intercept and two cadence indicators.\"\"\"\n", " n = len(data)\n", " return np.column_stack(\n", " [\n", " np.ones(n),\n", " np.full(n, float(arm == \"medium\")),\n", " np.full(n, float(arm == \"high\")),\n", " ]\n", " )\n", "\n", "\n", "saturated = MSM(\n", " design=cadence_indicators,\n", " terms=(\"(intercept)\", \"medium vs low\", \"high vs low\"),\n", " design_kind=\"known\",\n", ")\n", "saturated_result = study.identify(MSMProjection(saturated)).estimate(method=method)\n", "arm_contrasts = study.identify(ATE(reference=\"low\")).estimate(method=method)\n", "\n", "pairs = (\n", " (\"msm[(intercept)]\", arm_result, \"ey[low]\"),\n", " (\"msm[medium vs low]\", arm_contrasts, \"ate[medium vs low]\"),\n", " (\"msm[high vs low]\", arm_contrasts, \"ate[high vs low]\"),\n", ")\n", "rows = []\n", "for name, counterpart, key in pairs:\n", " coefficient, match = saturated_result[name], counterpart[key]\n", " rows.append(\n", " {\n", " \"coefficient\": name,\n", " \"psi\": coefficient.psi,\n", " \"ci\": f\"({coefficient.ci[0]:.4f}, {coefficient.ci[1]:.4f})\",\n", " \"counterpart\": key,\n", " \"counterpart psi\": match.psi,\n", " \"counterpart ci\": f\"({match.ci[0]:.4f}, {match.ci[1]:.4f})\",\n", " \"influence curves agree within 1e-12\": bool(\n", " np.allclose(coefficient.influence_curve, match.influence_curve, rtol=0, atol=1e-12)\n", " ),\n", " }\n", " )\n", "print(pd.DataFrame(rows).round(4).to_string(index=False))" ] }, { "cell_type": "markdown", "id": "saturated-reading", "metadata": {}, "source": [ "**What this output tells you.** The intercept is -0.0038, the same value as the `low` mean in\n", "Step 6. The other coefficients are 0.6437 and 1.4681, the same as the `medium` versus `low` and\n", "`high` versus `low` contrasts. Each pair of influence curves agrees within `1e-12` at every row, so\n", "the intervals match too.\n", "\n", "A saturated projection is a reparameterization of the per-cadence report, not a different analysis.\n", "It is also the fallback for a board that does not want to commit to a contact spacing.\n", "\n", "A saturated check cannot see a component that vanishes when the model fits. The projection weight\n", "and a link's curvature term are two examples. The\n", "[technical entry](../technical-reference/msm-projections.md#validation-issues-special-to-this-method)\n", "describes the checks that do test them." ] }, { "cell_type": "markdown", "id": "assessment-heading", "metadata": {}, "source": [ "## Step 12: diagnostics, what the fit can check\n", "\n", "`include_retargets=True` runs the truncation curve with the assessment. Truncation clips each fitted\n", "cadence probability into a bound, so no row gets an extreme weight. The curve retargets each\n", "coefficient at a range of bounds without refitting the nuisance models. The support report describes\n", "positivity in the fitted treatment model. See\n", "[diagnostics](../user-guide/results-assessment.md#diagnostics)." ] }, { "cell_type": "code", "execution_count": 12, "id": "assessment", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:24.032011Z", "iopub.status.busy": "2026-10-16T01:17:24.031718Z", "iopub.status.idle": "2026-10-16T01:17:24.210243Z", "shell.execute_reply": "2026-10-16T01:17:24.209113Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " surface check status\n", " validation score_equations passed\n", " validation support warning\n", " validation nuisance_models completed\n", "diagnostics truncation_curve completed\n", "diagnostics refute deferred\n", "needs attention: ('support',)\n", "\n", "Positivity / overlap diagnostics\n", "--------------------------------\n", "n = 3000; propensity truncated to [0.0114, 0.9886]\n", "\n", "g(W) quantile 0% 1% 5% 25% 50% 75% 95% 99% 100% \n", "------------- ------ ------ ------ ------ ------ ------ ------ ------ ------\n", "g[high] 0.0006 0.0153 0.0478 0.1684 0.3291 0.5263 0.7641 0.8681 0.9546\n", "g[low] 0.0148 0.0799 0.1332 0.2191 0.2755 0.3260 0.4022 0.4494 0.5720\n", "g[medium] 0.0044 0.0168 0.0470 0.1703 0.3416 0.5430 0.7773 0.8897 0.9847\n", "\n", "arm n effective n ESS / n top 1% weight top 5% weight\n", "------ ---- ----------- ------- ------------- -------------\n", "high 1079 514.9 0.477 0.086 0.218 \n", "low 816 727.2 0.891 0.032 0.105 \n", "medium 1105 655.3 0.593 0.068 0.190 \n", "\n", "threshold P(g < t) P(g > 1-t)\n", "--------- -------- ----------\n", "0.01 0.0033 0.0000 \n", "0.025 0.0136 0.0001 \n", "0.05 0.0371 0.0009 \n", "0.1 0.0963 0.0041 \n", "\n", "truncated: 36 unit(s) (1.20%); most extreme untruncated g(W) = 0.011537\n", "\n", "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 \n", "----- -------- ----------- -------------------- ------------- ------------ ----------- ----------- --------- --------- ---------------- ----------\n", "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%)\n", "(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)\n", "max |clever covariate| (msm): 268.5\n", "\n", "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.\n", "\n", "Score-equation check\n", "--------------------\n", "target kind |score| before threshold ratio ok \n", "---------------------- --------------- --------- --------- --------- -------- ---\n", "msm fluctuation 0.000e+00 5.344e-03 6.760e-07 0.00e+00 yes\n", "msm[(intercept)] influence curve 0.000e+00 - 6.760e-07 0.00e+00 yes\n", "msm[assigned contacts] influence curve 0.000e+00 - 1.569e-07 0.00e+00 yes\n", "\n", "PASS: the targeting step solved the estimated efficient score equation.\n", "\n", " bound psi delta_from_fitted ci_lower ci_upper truncated_fraction\n", "0.0010 0.2690 0.0002 0.2521 0.2858 0.0003\n", "0.0050 0.2689 0.0001 0.2520 0.2857 0.0023\n", "0.0100 0.2688 0.0000 0.2519 0.2856 0.0100\n", "0.0114 0.2688 0.0000 0.2519 0.2856 0.0120\n", "0.0250 0.2687 -0.0001 0.2519 0.2855 0.0403\n", "0.0500 0.2700 0.0013 0.2536 0.2865 0.1080\n", "0.1000 0.2699 0.0012 0.2546 0.2852 0.2677\n", "0.2000 0.2654 -0.0033 0.2525 0.2784 0.6093\n", "largest absolute slope movement: 0.0033\n" ] } ], "source": [ "assessment = trend_result.assess(include_retargets=True)\n", "ledger = assessment.to_frame()\n", "ran = ledger[~ledger[\"status\"].isin([\"unavailable\", \"not_applicable\"])]\n", "print(ran[[\"surface\", \"check\", \"status\"]].to_string(index=False))\n", "print(\"needs attention:\", tuple(item.name for item in assessment.attention))\n", "\n", "support = assessment.report(\"support\")\n", "scores = assessment.report(\"score_equations\")\n", "curve = assessment.report(\"truncation_curve\")\n", "print()\n", "print(support.summary())\n", "print()\n", "print(scores.summary())\n", "print()\n", "slope_curve = curve.loc[\n", " curve[\"estimand\"] == \"msm[assigned contacts]\",\n", " [\"bound\", \"psi\", \"delta_from_fitted\", \"ci_lower\", \"ci_upper\", \"truncated_fraction\"],\n", "]\n", "print(slope_curve.round(4).to_string(index=False))\n", "print(f\"largest absolute slope movement: {slope_curve['delta_from_fitted'].abs().max():.4f}\")" ] }, { "cell_type": "markdown", "id": "assessment-reading", "metadata": {}, "source": [ "**What this output tells you.** Read the parts in order.\n", "\n", "| output part | what it shows on this draw |\n", "| --- | --- |\n", "| 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 |\n", "| 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 |\n", "| 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 |\n", "| score-equation check | one row for the fluctuation's two-coefficient score, and one influence-curve row per coefficient |\n", "| truncation curve | the slope stays between 0.2654 and 0.2700 across bounds from 0.0010 to 0.2000 |\n", "\n", "The largest bound clips a fraction of 0.6093 of the rows, and the slope moves by at most 0.0033. On\n", "this draw, the slope depends little on the truncation bound. That does not verify positivity.\n", "\n", "The support report describes the fitted treatment model. It cannot show how small the true\n", "probabilities are. About one patient in a hundred has a true probability of `high` or\n", "`medium` below 0.01. The probe script `truth.py` gives the true probabilities." ] }, { "cell_type": "markdown", "id": "sensitivity-heading", "metadata": {}, "source": [ "## Step 13: sensitivity on the cadence contrasts\n", "\n", "Sensitivity analysis asks how strong an unmeasured confounder would need to be to change the\n", "conclusion. [Sensitivity analysis](../user-guide/results-assessment.md#sensitivity-analysis) defines\n", "it. `cleverly` implements the omitted-variable bound for arm-indexed estimands, not for an MSM\n", "coefficient. The code reads the robustness value of each contrast of the `ATE(reference=\"low\")` fit\n", "from Step 11. Those contrasts have the same influence curves as the saturated coefficients." ] }, { "cell_type": "code", "execution_count": 13, "id": "sensitivity", "metadata": { "execution": { "iopub.execute_input": "2026-10-16T01:17:24.212453Z", "iopub.status.busy": "2026-10-16T01:17:24.211761Z", "iopub.status.idle": "2026-10-16T01:17:24.219671Z", "shell.execute_reply": "2026-10-16T01:17:24.218805Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ate[medium vs low]: robustness value 0.204 (confidence-limit value 0.181)\n", "ate[high vs low]: robustness value 0.420 (confidence-limit value 0.394)\n" ] } ], "source": [ "robustness = {\n", " name: arm_contrasts.sensitivity.robustness_value(estimand=name)\n", " for name in (\"ate[medium vs low]\", \"ate[high vs low]\")\n", "}\n", "for name, values in robustness.items():\n", " print(\n", " f\"{name}: robustness value {values['rv']:.3f} (confidence-limit value {values['rva']:.3f})\"\n", " )" ] }, { "cell_type": "markdown", "id": "sensitivity-reading", "metadata": {}, "source": [ "**What this output tells you.** Each robustness value assumes worst-case alignment (`rho=1`). It\n", "is the equal outcome-side and treatment-side strength that moves the point estimate to zero. The\n", "confidence-limit value is the strength at which the confidence limit reaches zero. Both strengths\n", "are partial R² values between 0 and 1.\n", "\n", "| contrast | robustness value | confidence-limit value |\n", "| --- | --- | --- |\n", "| `ate[medium vs low]` | 0.204 | 0.181 |\n", "| `ate[high vs low]` | 0.420 | 0.394 |\n", "\n", "The `medium` versus `low` contrast reaches zero at the weaker confounding strength. These values\n", "bound the cadence contrasts, not the slope. The slope is a fixed contrast of the three means\n", "(Step 7), but this bound does not cover that combination. The\n", "[omitted-variable bounds](../technical-reference/validation-methods.md#omitted-variable-bounds-robustness-value-benchmark-and-contours)\n", "section defines each quantity.\n", "\n", "The bound needs a finite second moment of the Riesz representer, so $E[1/g]$ must be finite for\n", "each cadence. On this law it is finite for all three. The review probe script `truth.py` computes it\n", "in closed form from the cadence logits." ] }, { "cell_type": "markdown", "id": "method-evidence", "metadata": {}, "source": [ "## Method evidence\n", "\n", "The fits on this page use one draw of a synthetic law. The registered studies in the table\n", "measure each method over repeated draws. Each study page gives its law, its nuisance models,\n", "and the result of each cell.\n", "\n", "| fit on this page | registered study | law of the study |\n", "| --- | --- | --- |\n", "| Step 6, each cadence's mean | [Ordinary multi-arm TMLE](../technical-reference/method-evidence/ordinary-multi-arm-tmle.md) | another law |\n", "| Step 7, the trend | [Point-treatment MSM projection](../technical-reference/method-evidence/point-treatment-msm-projection.md) | another law |\n", "\n", "No registered study covers three cadences, fitted nuisances, or simultaneous bands for a\n", "projection.\n", "\n", "A review probe is a script under `reviews/notebook-review/probes/`. The step readings quote its\n", "results. A review probe is not a registered study.\n", "\n", "The [technical entry](../technical-reference/msm-projections.md#validation-issues-special-to-this-method) lists the evidence and the limits of the method." ] }, { "cell_type": "markdown", "id": "interpret", "metadata": {}, "source": [ "## What the analysis tells the board\n", "\n", "| question | answer from this analysis |\n", "| --- | --- |\n", "| 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 |\n", "| the population and the intervention | adults with a discharge-home order, each assigned one cadence and analyzed as assigned |\n", "| 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 |\n", "| whether one line is adequate | no. The miss at `medium` has the interval (-0.2566, -0.1602), which excludes zero |\n", "| what supports the reading | the score-equation check, the saturated control, and the stable truncation curve |\n", "| 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 |\n", "| what supports the cadence contrasts | their robustness values, 0.204 and 0.420. These values do not cover the slope |\n", "| 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 |\n", "\n", "Report the three cadence means or the two contrasts beside the slope. The gain per contact is not\n", "constant, and the slope hides that. The analysis does not support these conclusions:\n", "\n", "| conclusion | why the analysis does not support it |\n", "| --- | --- |\n", "| 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 |\n", "| ten contacts would raise the score further | no patient was assigned more than six contacts |\n", "| the board should choose the `high` cadence | that choice needs the cadence contrasts, their costs, and the board's own trade-offs |\n", "| the slope is robust to unmeasured confounding | the omitted-variable bound covers the cadence contrasts, not the slope |\n", "| the interval covers at 95% in general | one review probe on one law measured 0.9405, and no registered study covers this configuration |" ] }, { "cell_type": "markdown", "id": "where-next", "metadata": {}, "source": [ "## Where to go next\n", "\n", "The same projection works over regimens and horizons in a longitudinal fit. Its design callable\n", "also receives the horizon, and it codes each regimen as a number in the same way as Step 8. A\n", "longitudinal MSM projection cross-fits by default, with one pooled update per node.\n", "[Longitudinal MSM projections](../user-guide/longitudinal.md#longitudinal-msm-projections) shows the\n", "call. [Sequential offers](sequential-offers.ipynb) and\n", "[time-to-event outcomes](longitudinal-survival.ipynb) estimate the per-regimen parameters that the\n", "projection summarizes.\n", "\n", "The [examples index](index.md#the-program) lists every tutorial in the program." ] } ], "metadata": { "cleverly_execution": { "command": [ "python", "scripts/execute_notebook.py", "docs/examples/msm-projections.ipynb" ], "gated": { "code_source_sha256": "65602e89402ddf0f5068c69eb59c9c415725db96b31efcbb8a0364e56a3bc4fc", "execution_payload_sha256": "e17fd4030308cc087cdf95da52a1dd0850e1b83a7b25408abd910b1c62665924" }, "recorded": { "cleverly_commit": "4cea0d1ad32abd00580cd9f3a27d84a8c294db5f", "cleverly_version": "0.1.3", "cleverly_worktree_clean": false, "dependency_lock_sha256": "2e4f7ad3264ff9b35cf7588e5c93a0e20fc5c0852092ac1da365491543836682", "generator_files": [ "scripts/execute_notebook.py", "tests/notebooks.py" ], "generator_sha256": "f247771ab390b23e7e092d7eb35794e241c2ac84fc0526695df1cee1fc457a83", "package_source_sha256": "5fce085883c189351f5a0d142b216264208c189d01a46bb372a5e8fc68eed390" }, "schema_version": 5 }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.16" } }, "nbformat": 4, "nbformat_minor": 5 }