23  Observational Methods: Learning from Data

Status: Draft

v0.4

23.1 Introduction

NotePart III — Observable

Stratum. What the observation process records, and how we estimate interventional and counterfactual targets under partial observation.

Levels. L1 association and validation; L2 G-methods / TMLE / LMTP / forecasting under \(do(\cdot)\); L3 unit-level counterfactuals, then hypothesis generation and redesign of measurement (Ch. 25–27).

Deferred. Full CDM notation and package end-to-end patterns (Part IV). Mechanistic ODE/\(do(\cdot)\) simulation already lived in Part II—here the same L2 questions are answered from observational panels.

How to read. Centripetal return: climb from the recorded surface back towards structural and dynamical hypotheses with explicit estimands, assumptions, certificates, and estimators. A fitted result does not by itself establish identifiability.

Recorded data can constrain causal structure and effects, but only through an explicit observation process. This chapter develops Seeing in the Observable role: using observable data to learn about structural and dynamical mechanisms while keeping the record-generating process separate from the process being measured.

23.2 Time-Varying Confounding

23.2.1 Definition

Time-varying confounding occurs when (Robins 1986; Robins et al. 2000):

  1. Past confounders \(L_t\) affect both treatment \(A_t\) and outcome \(Y_t\)
  2. Past treatment \(A_{t-1}\) affects confounders \(L_t\)
  3. Confounders \(L_t\) affect future treatment \(A_{t+1}\)

Result: Standard methods (e.g., regression adjustment) fail because confounding evolves over time.

23.2.2 Example: Treatment and Disease Severity

  • Treatment \(A_t\): Drug at time \(t\)
  • Confounder \(L_t\): Disease severity at time \(t\)
  • Outcome \(Y_t\): Recovery at time \(t\)

Problem:

  • \(L_t\) affects \(A_t\) (severity influences treatment choice)
  • \(L_t\) affects \(Y_t\) (severity influences recovery)
  • \(A_{t-1}\) affects \(L_t\) (treatment affects severity)

This creates a feedback loop: treatment affects confounder, which affects future treatment.

23.2.3 Framing in CDMs

A CDM naturally handles time-varying confounding:

\[ \begin{aligned} L_t &\coloneqq f_L(L_{t-1}, A_{t-1}, Y_{t-1}, U^L_t) \quad \text{(confounder dynamics)} \\ A_t &\coloneqq \pi(H_t, U^A_t) \quad \text{(treatment assignment)} \\ Y_t &\coloneqq f_Y(L_t, A_t, Y_{t-1}, U^Y_t) \quad \text{(outcome)} \end{aligned} \]

The CDM structure makes the confounding explicit and allows us to adjust for it properly.

23.3 G-Methods: Handling Time-Varying Confounding

G-methods are a class of methods designed for time-varying treatments and confounders (Robins 1986; Robins et al. 2000; Hernán and Robins 2020). The “G” stands for “generalised” or “g-computation.”

23.3.1 G-Computation (G-Formula)

Idea: Simulate what would happen under a treatment strategy by integrating over confounder distributions.

Formula: \[ \mathbb{E}[Y^{do(\bar{a}_{1:T})}] = \int \mathbb{E}[Y \mid \bar{a}_{1:T}, \bar{l}_{1:T}] \prod_{t=1}^{T} P(l_t \mid \bar{a}_{t-1}, \bar{l}_{t-1}) \, d\bar{l}_{1:T} \]

Steps:

  1. Fit models for confounder evolution: \(P(L_t \mid \bar{A}_{t-1}, \bar{L}_{t-1})\)
  2. Fit outcome model: \(\mathbb{E}[Y \mid \bar{A}_{1:T}, \bar{L}_{1:T}]\)
  3. Simulate confounder trajectories under treatment strategy
  4. Compute expected outcomes

Mathematical form: The g-formula expresses the interventional distribution as:

\[ P^{do(\bar{A}_t = \bar{a}_t)}(Y_t) = \int \prod_{s=1}^t P(Y_s \mid \bar{A}_s = \bar{a}_s, \bar{L}_s = \bar{l}_s) \prod_{s=1}^t P(L_s \mid \bar{A}_{s-1} = \bar{a}_{s-1}, \bar{L}_{s-1} = \bar{l}_{s-1}) \, d\bar{l}_t \]

Advantages:

  • Handles time-varying confounding naturally
  • Can estimate effects of complex treatment strategies
  • Provides full distribution, not just means

Limitations:

  • Requires correct specification of outcome and confounder models
  • Sensitive to model misspecification
  • Computationally intensive for long time horizons

23.3.2 Implementation: G-Computation

The example in the chunk-g-computation chunk fits both transitions to an observational cohort. At each step \(t\), treatment \(A_{t-1}\) is chosen from the previous severity score \(L_{t-1}\); severity \(L_t\) and outcome \(Y_t\) then follow. The continuous, centred severity score is not restricted to \([0,1]\). Identification assumes consistency, sequential exchangeability given the recorded history, and positive probability of both treatments. The linear Gaussian transition models are correctly specified for this simulation (Robins 1986).

The shared implementation fits coefficients by least squares. It samples initial states from the observed cohort and simulates always-treat and never-treat policies. The difference in means estimates the terminal-outcome contrast; the outcome histograms show individual variation, not confidence intervals for that contrast. Chapter 28 adds uncertainty and model criticism to a related fitted example.

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using StableRNGs CairoMakie
include(joinpath(dirname(Base.active_project()), "scripts", "longitudinal_gcomputation.jl"))

T = 5
data = simulate_treatment_cohort(1000; T=T, rng=StableRNG(42))
confounder_model = fit_confounder_model(data)
outcome_model = fit_outcome_model(data)
baseline = [(s.L[1], s.Y[1]) for s in data]
strategy_always_treat = (t, L, Y) -> 1.0
strategy_never_treat = (t, L, Y) -> 0.0
# Common random numbers reduce Monte Carlo variance of the policy comparison.
outcomes_treat = g_compute(strategy_always_treat, confounder_model, outcome_model,
    5000; T=T, baseline=baseline, rng=StableRNG(123))
outcomes_no_treat = g_compute(strategy_never_treat, confounder_model, outcome_model,
    5000; T=T, baseline=baseline, rng=StableRNG(123))
println("Fitted confounder transition: ", confounder_model)
println("Fitted outcome transition: ", outcome_model)
println("Terminal mean, always treat: ", round(mean(outcomes_treat); digits=3))
println("Terminal mean, never treat: ", round(mean(outcomes_no_treat); digits=3))
println("Estimated policy contrast: ", round(mean(outcomes_treat - outcomes_no_treat); digits=3))
Fitted confounder transition: (α = 0.09903671285464175, β = -0.20196445144385164, γ = 0.795981953512425, σ = 0.09958645169367795)
Fitted outcome transition: (α = -0.0037585903543576526, β = 0.5022426750522038, γ = -0.2901796412286034, δ = 0.6995986330011287, σ = 0.09961561961206072)
Terminal mean, always treat: 1.247
Terminal mean, never treat: -0.368
Estimated policy contrast: 1.615
Figure 23.1: G-computation: simulating outcomes under treatment strategy

23.3.3 Inverse Probability of Treatment Weighting (IPTW)

Idea: Weight observations by inverse probability of receiving their observed treatment, creating a pseudo-population where treatment is independent of confounders.

Weights: \[ w_i = \prod_{t=1}^{T} \frac{1}{P(A_t = a_{it} \mid \bar{A}_{t-1}, \bar{L}_{t})} \]

Use: Weighted regression or weighted means to estimate causal effects.

Stabilised weights (reduce variance): \[ w_i^{\text{stab}} = \prod_{t=1}^{T} \frac{P(A_t = a_{it} \mid \bar{A}_{t-1})}{P(A_t = a_{it} \mid \bar{A}_{t-1}, \bar{L}_{t})} \]

Interpretation: Subjects with treatment histories that are unlikely given their confounder history receive higher weights, creating a weighted population where treatment is independent of confounders.

Steps:

  1. Fit treatment models: Estimate \(P(A_s \mid \bar{A}_{s-1}, \bar{L}_s)\) for each time \(s\)
  2. Compute weights: Calculate IPTW weights \(w_t\) for each subject
  3. Stabilize weights (optional): Use stabilized weights to reduce variance
  4. Fit outcome model: Fit a model for \(E[Y_t \mid \bar{A}_t]\) on weighted data

Advantages:

  • Intuitive: creates a pseudo-randomized population
  • Can use standard regression on weighted data
  • Handles time-varying confounding

Limitations:

  • Requires correct specification of treatment model
  • Extreme weights can create instability
  • Sensitive to positivity violations

23.3.4 Implementation: IPTW

The following example implements IPTW for time-varying treatments:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using Random Distributions GLM DataFrames CairoMakie

Random.seed!(42)

# Simulate data with time-varying confounding
n_subjects = 200
T = 5
data = []

for i in 1:n_subjects
    L = [rand(Normal(0.5, 0.2))]  # Initial severity
    A = Int[]
    Y = [rand(Normal(0.0, 0.1))]  # Initial outcome
    
    for t in 2:T
        # Treatment assignment (depends on L_{t-1} - confounding!)
        p_treat = 1 / (1 + exp(-2*(L[end] - 0.5)))
        push!(A, rand(Bernoulli(p_treat)))
        
        # Confounder evolution
        L_new = 0.1 - 0.2*A[end] + 0.8*L[end] + rand(Normal(0, 0.1))
        push!(L, max(0, min(1, L_new)))
        
        # Outcome
        Y_new = 0.5*A[end] - 0.3*L[end] + 0.7*Y[end] + rand(Normal(0, 0.1))
        push!(Y, Y_new)
    end
    
    push!(data, (L=L, A=A, Y=Y))
end

# Step 1: Fit treatment models P(A_t | A_{t-1}, L_t)
# For each time point, fit logistic regression
treatment_models = []

for t in 2:T
    # Prepare data for time t
    A_t = [d.A[t-1] for d in data]  # A_t
    A_prev = [t > 2 ? d.A[t-2] : 0 for d in data]  # A_{t-1}
    L_t = [d.L[t-1] for d in data]  # L_t
    
    df = DataFrame(A_t = A_t, A_prev = A_prev, L_t = L_t)
    
    # Fit logistic regression
    # NOTE: In production, use flexible models (splines, machine learning) for robustness
    # P(A_t = 1 | A_{t-1}, L_t) = logistic(α + β*A_{t-1} + γ*L_t)
    model = glm(@formula(A_t ~ A_prev + L_t), df, Binomial(), LogitLink())
    push!(treatment_models, model)
end

# Step 2: Compute IPTW weights
# w_i = ∏_t 1 / P(A_t = a_{it} | A_{t-1}, L_t)
weights = Float64[]

for d in data
    w = 1.0
    for t in 2:T
        A_t = d.A[t-1]
        A_prev = t > 2 ? d.A[t-2] : 0
        L_t = d.L[t-1]
        
        # Predict P(A_t = a_{it} | A_{t-1}, L_t)
        model = treatment_models[t-1]
        df_t = DataFrame(A_prev = [A_prev], L_t = [L_t])
        p_pred = GLM.predict(model, df_t)[1]
        
        # Weight = 1 / probability of observed treatment
        # Protect against division by zero
        p_observed = A_t == 1 ? max(p_pred, 1e-10) : max(1 - p_pred, 1e-10)
        w *= 1.0 / p_observed
    end
    push!(weights, w)
end

# Step 3: Stabilise weights (optional but recommended)
# w_stab = ∏_t P(A_t | A_{t-1}) / P(A_t | A_{t-1}, L_t)
# Numerator: marginal treatment probability (no confounders)
stabilised_weights = Float64[]

for d in data
    w_stab = 1.0
    for t in 2:T
        A_t = d.A[t-1]
        A_prev = t > 2 ? d.A[t-2] : 0
        L_t = d.L[t-1]
        
        # Marginal model: P(A_t | A_{t-1}) only
        # Fit on all data at this time point (not just this subject)
        A_t_all = [d.A[t-1] for d in data]
        A_prev_all = [t > 2 ? d.A[t-2] : 0 for d in data]
        df_marginal_all = DataFrame(A_t = A_t_all, A_prev = A_prev_all)
        model_marginal = glm(@formula(A_t ~ A_prev), df_marginal_all, Binomial(), LogitLink())
        p_marginal = GLM.predict(model_marginal, DataFrame(A_prev = [A_prev]))[1]
        
        # Conditional model: P(A_t | A_{t-1}, L_t)
        model_cond = treatment_models[t-1]
        df_cond = DataFrame(A_prev = [A_prev], L_t = [L_t])
        p_cond = GLM.predict(model_cond, df_cond)[1]
        
        # Stabilised weight
        # Protect against division by zero
        p_observed = A_t == 1 ? max(p_cond, 1e-10) : max(1 - p_cond, 1e-10)
        p_marginal_observed = A_t == 1 ? max(p_marginal, 1e-10) : max(1 - p_marginal, 1e-10)
        w_stab *= p_marginal_observed / p_observed
    end
    push!(stabilised_weights, w_stab)
end

# Step 4: Estimate treatment effect using weighted outcome model
# Compare outcomes under always treat vs never treat
final_outcomes = [d.Y[end] for d in data]
always_treat = [all(a .== 1) for a in [d.A for d in data]]

# Weighted mean difference
sum_w_treat = sum(stabilised_weights .* always_treat)
sum_w_no_treat = sum(stabilised_weights .* (.!always_treat))
mean_treat = sum_w_treat > 0 ? sum(stabilised_weights .* final_outcomes .* always_treat) / sum_w_treat : 0.0
mean_no_treat = sum_w_no_treat > 0 ? sum(stabilised_weights .* final_outcomes .* (.!always_treat)) / sum_w_no_treat : 0.0

println("IPTW results:")
println("  Mean unstabilised weight: ", round(mean(weights), digits=2))
println("  Mean stabilised weight: ", round(mean(stabilised_weights), digits=2))
println("  Max unstabilised weight: ", round(maximum(weights), digits=2))
println("  Max stabilised weight: ", round(maximum(stabilised_weights), digits=2))
println("\nTreatment effect (weighted):")
println("  E[Y | always treat] = ", round(mean_treat, digits=3))
println("  E[Y | never treat] = ", round(mean_no_treat, digits=3))
println("  Effect = ", round(mean_treat - mean_no_treat, digits=3))
IPTW results:
  Mean unstabilised weight: 16.71
  Mean stabilised weight: 1.03
  Max unstabilised weight: 147.51
  Max stabilised weight: 9.19

Treatment effect (weighted):
  E[Y | always treat] = 1.166
  E[Y | never treat] = 0.308
  Effect = 0.859
Figure 23.2: IPTW: creating pseudo-population with inverse probability weights

23.3.5 Doubly Robust Methods

Idea: Combine outcome model and treatment model. Consistent if either model is correct.

Examples:

23.3.6 Variational backdoor head (RxInfer)

When the estimand is a scalar total effect and confounders are known, CausalDynamics can residualise and hand a conjugate GraphPPL model to RxInfer (see Chapter 5 and Chapter 30):

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using CausalDynamics Graphs DataFrames Random
using RxInfer

Random.seed!(7)
n = 800
Z = randn(n)
X = 0.8 .* Z .+ randn(n)
Y = 1.5 .* X .+ 0.6 .* Z .+ 0.4 .* randn(n)
g = DiGraph(3)
add_edge!(g, 1, 2); add_edge!(g, 1, 3); add_edge!(g, 2, 3)
df = DataFrame(Z = Z, X = X, Y = Y)
t_o = @elapsed begin
    global ate_rx = infer_backdoor_effect(g, df, 2, 3;
        node_names = Dict(1 => :Z, 2 => :X, 3 => :Y), iterations = 30)
end
println("RxInfer ATE τ_mean ≈ ", round(ate_rx.τ_mean; digits = 3),
    " (truth 1.5);  ", round(t_o; digits = 3), " s")
RxInfer ATE τ_mean ≈ 1.475 (truth 1.5);  8.515 s

Use this when you want a fast posterior on \(\tau\) under a linear Gaussian head. Keep TMLE / CausalTargeted for doubly robust or continuous MTP estimands.

23.4 Challenges in Time-Varying Settings

23.4.1 Positivity

Positivity (also called “overlap” or “common support”) requires that all treatment histories have positive probability given confounder history:

\[ P(A_t = a_t \mid \bar{A}_{t-1} = \bar{a}_{t-1}, \bar{L}_t = \bar{l}_t) > 0 \]

for all possible treatment values \(a_t\) and all confounder histories \(\bar{l}_t\) that occur in the data.

Positivity violation: Some treatment histories have zero (or near-zero) probability given confounder history.

Problem:

  • Cannot estimate effects for impossible treatment histories
  • IPTW weights become infinite (or very large)
  • G-computation cannot simulate under impossible interventions
  • Estimates become unstable or undefined

Solutions:

  1. Restrict to feasible policies: Only estimate effects for treatment strategies that are possible given the data
  2. Truncate weights: Cap IPTW weights at a maximum value (e.g., 10 or 20) to reduce instability
  3. Use bounds: Compute bounds on treatment effects when positivity is violated
  4. Sensitivity analysis: Assess how results change with different positivity assumptions
  5. Continuous MTP diagnostics: For additive shifts of continuous exposures, CausalTargeted.jl positivity_report (and grid positivity=true) summarises support / clamp status before interpreting TE curves (Petersen et al. 2012; Díaz et al. 2023), see Policy Evaluation

23.4.2 Model Misspecification

Problem: Models for confounders, treatment, or outcome may be incorrectly specified, leading to biased estimates.

G-computation: Sensitive to misspecification of outcome and confounder models. If these models are wrong, estimates will be biased.

MSMs (IPTW): Sensitive to misspecification of treatment model. If the treatment model is wrong, weights will be incorrect and estimates biased. For repeated outcomes under nested sampling units (farms, classrooms), CausalTargeted’s repeated-outcome MSM can keep the same IF point estimates and pass cluster= so \(\widehat{\Sigma}\) uses a cluster sandwich (Ch. 28; Ch. 30).

Solutions:

  1. Use flexible models: Machine learning methods (random forests, neural networks, Super Learner (Laan et al. 2007)) can reduce misspecification
  2. Model criticism: Iterate simulate-and-recover, calibration, residuals, and posterior predictive checks (see Model Validation with Observable Data; (Gelman et al. 2026))
  3. Sensitivity analysis: Assess how results change with different model specifications
  4. Robust methods: TMLE provides some protection through double robustness (see TMLE and Doubly Robust Estimation)

23.4.3 Unmeasured Confounding

Problem: Unobserved confounders bias estimates even when all measured confounders are adjusted for.

Time-varying setting: Unmeasured time-varying confounders are particularly problematic because they can create complex dependencies that are difficult to detect.

Solutions:

  1. Sensitivity analysis: Quantify how robust results are to unmeasured confounding (see Identification: When Can We Learn from Data? for E-values and sensitivity parameters; CausalTargeted sensitivity_report for continuous MTP tipping-point / partial-\(R^2\) diagnostics (Cinelli and Hazlett 2020))
  2. Bounds: Compute worst-case bounds on treatment effects under assumptions about unmeasured confounders
  3. Instrumental variables: Use instruments that are independent of unmeasured confounders (when available)
  4. Negative controls: Use known null effects to detect unmeasured confounding
  5. Design improvements: Collect data on suspected unmeasured confounders in future studies

23.5 Incomplete observation

What is recorded is never the whole of what became. An assay \(Y_t\) is a partial observation of latent state; when the assay itself is absent, or a covariate fails to land in the table, the Observable surface thins further. Incomplete observation belongs inside the measurement story: it is part of how a design records (or fails to record) a process of becoming, and it is organised on the same two grids that structure this book, strata (where the incompleteness lives) and levels of Reason (what query we are asking) (Pearl 2009; Laan and Rose 2011; Little and Rubin 2019).

Notation is fixed in the Concept Reference: complete \(V\), response \(R_V\), recorded \(V^{\mathrm{rec}}\) (figure alias \(V^*\)), fills \(\tilde{V}\) off the causal graph, counterfactuals \(V^{do(\cdot)}(\mathbf{u})\). The m-DAG atlas for MCAR–MNAR is in Chapter 18; the same stratum × level placement is argued here for estimation and for the Dynamical fork between censoring and MAR missing terminals.

23.5.1 The same gap, three strata

Structural. Response indicators are themselves mechanisms. For an outcome \(Y\):

\[ \begin{aligned} Y &\coloneqq f_Y\bigl(\mathrm{pa}_G(Y),\, U^Y\bigr), \\ R_Y &\coloneqq f_{R_Y}\bigl(\mathrm{pa}_G(R_Y),\, U^{R_Y}\bigr), \\ Y^{\mathrm{rec}} &\coloneqq \begin{cases} Y & R_Y = 1, \\ \texttt{missing} & R_Y = 0. \end{cases} \end{aligned} \]

MCAR, MAR, and MNAR name different Structural claims about the parents of \(R_Y\) and about \(U^{R_Y}\) on \(G\). CausalDynamics states that claim as a MissingnessSpec, certifies it (certify_missingness), and may attach the certificate when identify(...; missingness=) runs. Causal identification of \(P(Y \mid do(A))\) and identification under incomplete observation are reported separately: an MNAR claim can leave the missingness certificate unidentified while the backdoor still holds for the complete-data query. Keeping those certificates distinct leaves the freeze-frame of \(G\) explicit about the observation process as well as about treatment.

Figure 23.3 pairs the usual confounded dyad with \(R_Y\): under MAR, \(W\to R_Y\); under MNAR, \(Y\to R_Y\). The backdoor set \(\{W\}\) for \(A\to Y\) is unchanged; only the missingness claim differs.

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using DAGMakie CairoMakie Graphs

# Nodes: 1=W, 2=A, 3=Y, 4=R_Y, 5=Y* (≡ Y^rec)
g_mar_ry = SimpleDiGraph(5)
add_edge!(g_mar_ry, 1, 2); add_edge!(g_mar_ry, 1, 3); add_edge!(g_mar_ry, 2, 3)
add_edge!(g_mar_ry, 3, 5); add_edge!(g_mar_ry, 4, 5); add_edge!(g_mar_ry, 1, 4)

g_mnar_ry = SimpleDiGraph(5)
add_edge!(g_mnar_ry, 1, 2); add_edge!(g_mnar_ry, 1, 3); add_edge!(g_mnar_ry, 2, 3)
add_edge!(g_mnar_ry, 3, 5); add_edge!(g_mnar_ry, 4, 5); add_edge!(g_mnar_ry, 3, 4)

# A sits above the W–Y chord. MAR: R_Y under W, Y* under Y (rectangle).
# MNAR: R_Y under Y, Y* to the right of Y.
layout_ry_mar = Point2f[
    Point2f(0.0, 1.0), Point2f(1.2, 1.85), Point2f(2.4, 1.0),
    Point2f(0.0, 0.0), Point2f(2.4, 0.0),
]
layout_ry_mnar = Point2f[
    Point2f(0.0, 1.0), Point2f(1.2, 1.85), Point2f(2.4, 1.0),
    Point2f(2.4, 0.0), Point2f(3.6, 1.0),
]
const _lab_ry = ["W", "A", "Y", "R_Y", "Y*"]
const _col_o = RGBf(0.30, 0.55, 0.70)
const _col_l = RGBf(0.85, 0.65, 0.05)
const _col_r = RGBf(0.55, 0.55, 0.55)
Figure 23.3: Backdoor confounding with incomplete outcome. Left: MAR (\(W\to R_Y\)). Right: MNAR (\(Y\to R_Y\)). \(Y^*\equiv Y^{\mathrm{rec}}\); gold \(Y\) is complete; grey \(R_Y\) is recording. Causal backdoor adjustment for \(A\to Y\) and the missingness certificate are separate claims.

Dynamical. Gaps have a clock. A missing terminal \(Y\) on a two-time panel, an intermittent \(Y_t\) along a sequential path, and right-censoring in discrete-time survival leave similar holes in a table while answering different dynamical questions. Generative oracles (apply_missingness_mechanism, simulate_incomplete_panel) mask complete trajectories under an explicit mechanism; CDM solvers leave gaps marked rather than inventing continuous-looking fills. Survival censoring IPCW and MAR missingness of a terminal event indicator \(S_T\) therefore carry distinct assumptions and distinct APIs, much as association and intervention do on Pearl’s ladder.

Figure 23.4 contrasts a terminal response \(R_{S_T}\) (MAR missing event indicator) with a censoring node \(C_t\) on the survival path: same table hole, different object.

# Nodes: left panel MAR S_T — 1=W, 2=A, 3=S_T, 4=R_ST, 5=S*
# Right panel censoring — 1=W, 2=A, 3=S_T, 4=C_t
g_mar_st = SimpleDiGraph(5)
add_edge!(g_mar_st, 1, 2); add_edge!(g_mar_st, 1, 3); add_edge!(g_mar_st, 2, 3)
add_edge!(g_mar_st, 3, 5); add_edge!(g_mar_st, 4, 5); add_edge!(g_mar_st, 1, 4)

g_cens = SimpleDiGraph(4)
add_edge!(g_cens, 1, 2); add_edge!(g_cens, 1, 3); add_edge!(g_cens, 2, 3)
add_edge!(g_cens, 1, 4); add_edge!(g_cens, 4, 3)

layout_st = Point2f[
    Point2f(0.0, 1.0), Point2f(1.2, 1.85), Point2f(2.4, 1.0),
    Point2f(0.0, 0.0), Point2f(2.4, 0.0),
]
# Diamond: A above, C_t below the W–S_T chord.
layout_c = Point2f[
    Point2f(0.0, 1.0), Point2f(1.3, 1.85), Point2f(2.6, 1.0), Point2f(1.3, 0.15),
]
Figure 23.4: Dynamical fork. Left: MAR missing terminal event indicator (\(W\to R_{S_T}\), recorded \(S^*\)). Right: right-censoring \(C_t\) on the survival path. Each object has its own IPCW semantics.

Observable. Julia missing is the recorded token of \(R_V=0\): the surface mark that something was not grasped. Estimation policies live here, in CausalTargeted and CausalMediation: complete-case :drop, outcome IPCW :ipcw, covariate mean-imputation :impute, their combination :ipcw_impute, and (opt-in) posterior draws impute_posterior under a MAR certificate, pooled across LMTP grids by Rubin’s rule. Fills enter analysis as \(\tilde{Y}\) or \(Y^{\mathrm{imp}}\) under a documented policy; Structural encode (encode_to_panel) requires complete assays (or an explicit prior drop), so incomplete high-dimensional matrices stay incomplete until an Observable step has acted.

23.5.2 Levels of Reason under incomplete records

At L1, association under MCAR or MAR can often proceed with explicit drop or reweighting, while remaining clear that complete cases are themselves a selected process. At L2, interventional targets with partial \(Y\) need the missingness DAG (or an equivalent claim) alongside \(do(A)\): IPCW and pattern-mixture style identification act on \(R_Y\) and on the observed cells. At L3, unit-level counterfactuals ask what an explicitly declared cross-world coupling implies under another \(do(\cdot)\), commonly with an exogenous history held fixed. When \(Y\) (or a parent on the path) was unobserved, inference over that history and over the missing cell is more uncertain: \(\tilde{Y}\) under a MAR/MNAR claim may support estimation, while unit-level \(Y^{do(a)}(\mathbf{u})\) still requires a declared coupling and an explicit identification story. Shared-history replay with that claim, or an honest partial-identification statement, is the L3 path available here.

23.5.3 Recoverability before imputation

A causal effect may be identified from complete data yet fail to be recoverable from the records actually available. Recoverability asks whether the target is determined by the observed-data distribution under the declared missingness graph. The MCAR–MAR–MNAR classification alone does not answer this question. Some MNAR graphs permit recovery of particular targets; others leave observationally indistinguishable complete-data models with different target values (Mohan and Pearl 2018).

Graph-guided imputation therefore begins with the target and the recording mechanisms, rather than with a generic table-completion rule. For example, if baseline nutrition predicts both worm burden and attendance at follow-up, those relations belong in the observation model. Under an appropriate MAR claim, multiple imputation can use treatment, nutrition and other observed predictors to draw plausible missing outcomes. The imputation model must also accommodate interactions and non-linearities needed by the analysis model; a subsequent doubly robust estimator does not automatically repair incompatible imputations (Dashti et al. 2024). When the target is not recoverable, further measurements, additional identifying assumptions, bounds or sensitivity analysis are needed.

23.5.4 Multiple robustness depends on the missing variable

Missing outcomes, missing exposures and missing confounders yield different observed-data models. Zhang and colleagues develop a triply robust estimator for missing exposure information that is consistent when any two of three parametric models are correctly specified: treatment assignment, outcome regression and the exposure-recording mechanism. Their result also requires MAR, no unmeasured confounding and positivity for both treatment and recording (Zhang et al. 2016). The robustness guarantee belongs to that estimator and its model specification; counting three fitted nuisance models is insufficient to establish it. Chapter 20 describes the different robustness structure of a standard TMLE for missing outcomes.

23.5.5 Randomisation-based inference with missing outcomes

Randomisation supplies a known treatment-assignment distribution, but selective outcome recording can still invalidate a complete-case comparison. Heng and colleagues construct imputation-assisted randomisation tests by repeating the imputation procedure within each reassignment allowed by the experimental design. Under their missingness assumptions, this gives finite-population-exact control of type-I error for Fisher’s sharp null of no effect on any participant (Heng et al. 2026).

Their framework allows recording to depend on measured or unmeasured baseline characteristics and on the outcomes themselves. Treatment may affect recording through those outcomes, but a separate direct effect of treatment on recording falls outside the stated guarantee. Under the sharp null, the recording pattern therefore remains invariant across assignments. This permits some MNAR mechanisms without requiring a correctly specified imputation model; validity comes from the design and the missingness restriction. Testing this sharp null is distinct from estimating an average effect, and confidence regions require the additional assumptions specified for test inversion (Heng et al. 2026).

23.5.6 Scope of the stack

Package documentation keeps a strictly causal–statistical voice; the book cites that catalogue and leaves the heavy ledgers there. In outline:

Stratum Package surface Scope boundary
Structural MissingnessSpec, certificates on identify MAR requires a stated conditioning set when \(R_Y\) depends on covariates; MNAR stays unidentified without further assumptions
Dynamical Masked panels; sequential / survival handle_missing Censoring \(C_t\) and MAR \(R_{S_T}\) use distinct IPCW objects
Observable handle_missing, opt-in impute_posterior (\(\tilde{Y}\)) Complete assays (or explicit drop) before encode_to_panel; no silent Missing→Float64 in Dynamics

Heavy ledgers and recovery checks live in package stress notebooks (missingness grid; posterior; Deep SCM estimation). API layout: Software Patterns. Methods pages: CausalDynamics Missingness, CausalTargeted Missingness.

23.6 Stratum context

G-methods (G-computation and IPTW) use observable data to estimate causal effects under time-varying confounding. Incomplete observation follows the same placement: structural claims about \(R\), dynamical timing of gaps, and observation policies on missing.

23.7 Inverse Problems in Causal Dynamics

Inverse problems involve learning model parameters from observational data (Rackauckas 2026). This is central to causal inference: we observe outcomes and want to learn the mechanisms (parameters) that generated them.

23.7.1 Connection to Identification

Parameter identifiability is a prerequisite for causal identification (see Identification: When Can We Learn from Data?). If parameters are not identifiable, we cannot learn causal mechanisms from data.

Key Questions:

  1. Identifiability: Can parameters be uniquely determined from data?
  2. Estimation: How do we estimate parameters from finite data?
  3. Uncertainty: How uncertain are our parameter estimates?

23.7.2 Regularisation for Ill-Posed Problems

When inverse problems are ill-posed (multiple solutions or unstable), regularisation helps:

  • L2 regularisation: Penalise large parameter values
  • L1 regularisation: Encourage sparsity (few non-zero parameters)
  • Prior knowledge: Incorporate domain knowledge as priors

Example: Regularised parameter estimation

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using OrdinaryDiffEq StochasticDiffEq Optim Random LinearAlgebra

# SIR model
function sir!(du, u, p, t)
    S, I, R = u
    β, γ = p
    du[1] = -β * S * I
    du[2] = β * S * I - γ * I
    du[3] = γ * I
end

# Generate synthetic data
Random.seed!(42)
u0 = [0.99, 0.01, 0.0]
p_true = (0.3, 0.1)
tspan = (0.0, 50.0)
t_obs = 0.0:1.0:50.0

prob_true = ODEProblem(sir!, u0, tspan, p_true)
sol_true = solve(prob_true, Tsit5(), saveat = t_obs)
σ_obs = 0.01
observed_data = [max.(0.0, min.(1.0, u .+ σ_obs .* randn(3))) for u in sol_true.u]

# Unregularised cost function (data fit)
function cost(p)
    prob = ODEProblem(sir!, u0, tspan, p)
    sol = solve(prob, Tsit5(), saveat = t_obs)
    if sol.retcode != :Success
        return Inf
    end
    return sum([sum((sol.u[i] .- observed_data[i]).^2) for i in 1:length(sol.u)])
end

# Regularised cost function
function cost_regularised(p, λ)
    # Data fit term
    data_fit = cost(p)
    
    # Regularisation term (L2)
    reg = λ * sum(p.^2)
    
    return data_fit + reg
end

# Initial guess
p0 = [0.2, 0.15]  # β, γ

# Optimise without regularisation
result_unreg = optimize(cost, p0, LBFGS())
p_unreg = result_unreg.minimizer

# Optimise with regularisation (λ = 0.1)
result_reg = optimize(p -> cost_regularised(p, 0.1), p0, LBFGS())
p_reg = result_reg.minimizer

println("Regularised parameter estimation:")
println("  True parameters: β = $(p_true[1]), γ = $(p_true[2])")
println("  Unregularised: β = $(round(p_unreg[1], digits=3)), γ = $(round(p_unreg[2], digits=3))")
println("  Regularised (λ=0.1): β = $(round(p_reg[1], digits=3)), γ = $(round(p_reg[2], digits=3))")
println("  Regularisation reduces parameter magnitude: |p_reg| = $(round(norm(p_reg), digits=3)) vs |p_unreg| = $(round(norm(p_unreg), digits=3))")
Regularised parameter estimation:
  True parameters: β = 0.3, γ = 0.1
  Unregularised: β = 0.2, γ = 0.15
  Regularised (λ=0.1): β = 0.2, γ = 0.15
  Regularisation reduces parameter magnitude: |p_reg| = 0.25 vs |p_unreg| = 0.25

23.7.3 Multiple Shooting Methods

For long time horizons, multiple shooting breaks the problem into segments:

  1. Divide time span into segments: \([t_0, t_1], [t_1, t_2], \ldots, [t_{n-1}, t_n]\)
  2. Solve each segment independently
  3. Add continuity constraints: \(u_i(t_i) = u_{i+1}(t_i)\)

Advantages:

  • More robust than single shooting
  • Each segment can be solved in parallel
  • Better for long time horizons

Implementation:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using OrdinaryDiffEq StochasticDiffEq Optim Random CairoMakie

# SIR model
function sir!(du, u, p, t)
    S, I, R = u
    β, γ = p
    du[1] = -β * S * I
    du[2] = β * S * I - γ * I
    du[3] = γ * I
end

# Generate synthetic data
Random.seed!(42)
u0_true = [0.99, 0.01, 0.0]
p_true = (0.3, 0.1)
tspan = (0.0, 100.0)  # Long time horizon
t_obs = 0.0:2.0:100.0

prob_true = ODEProblem(sir!, u0_true, tspan, p_true)
sol_true = solve(prob_true, Tsit5(), saveat = t_obs)
σ_obs = 0.01
observed_data = [max.(0.0, min.(1.0, u .+ σ_obs .* randn(3))) for u in sol_true.u]

# Multiple shooting: divide into 4 segments
n_segments = 4
t_segments = [0.0, 25.0, 50.0, 75.0, 100.0]
segment_states = [Vector{Float64}(undef, 3) for _ in 1:(n_segments-1)]  # States at segment boundaries

# Cost function with continuity constraints
function cost_multiple_shooting(p_and_states)
    p = (p_and_states[1], p_and_states[2])
    total_cost = 0.0
    
    # Solve each segment
    for i in 1:n_segments
        t_start = t_segments[i]
        t_end = t_segments[i+1]
        
        # Initial condition for this segment
        if i == 1
            u0_seg = u0_true
        else
            u0_seg = [p_and_states[2 + 3*(i-2) + j] for j in 1:3]
        end
        
        prob_seg = ODEProblem(sir!, u0_seg, (t_start, t_end), p)
        sol_seg = solve(prob_seg, Tsit5(), saveat = t_obs[t_obs .>= t_start .&& t_obs .<= t_end])
        
        if sol_seg.retcode != :Success
            return Inf
        end
        
        # Data fit for this segment
        for (j, t) in enumerate(t_obs)
            if t_start <= t <= t_end
                idx = findfirst(==(t), sol_seg.t)
                if idx !== nothing && j <= length(observed_data)
                    total_cost += sum((sol_seg.u[idx] .- observed_data[j]).^2)
                end
            end
        end
        
        # Continuity constraint (except for last segment)
        if i < n_segments
            final_state = sol_seg.u[end]
            next_initial = [p_and_states[2 + 3*(i-1) + j] for j in 1:3]
            total_cost += 100.0 * sum((final_state .- next_initial).^2)  # Penalty for discontinuity
        end
    end
    
    return total_cost
end

# Initial guess: parameters + segment boundary states
p0 = [0.2, 0.15]  # β, γ
states0 = [u0_true for _ in 1:(n_segments-1)]
x0 = vcat(p0, vcat(states0...))

# Optimise
result = optimize(cost_multiple_shooting, x0, LBFGS())

# Extract estimated parameters
p_estimated = (result.minimizer[1], result.minimizer[2])

println("Multiple shooting results:")
println("  True parameters: β = $(p_true[1]), γ = $(p_true[2])")
println("  Estimated: β = $(round(p_estimated[1], digits=3)), γ = $(round(p_estimated[2], digits=3))")
Multiple shooting results:
  True parameters: β = 0.3, γ = 0.1
  Estimated: β = 0.2, γ = 0.15

When to Use Multiple Shooting:

  • Long time horizons where single shooting fails
  • Poor initial parameter guesses
  • Need for parallel computation

23.7.4 Collocation Methods

Collocation methods fit splines to data and match derivatives, avoiding ODE solving:

  1. Fit spline \(\tilde{u}(t)\) to observations
  2. Estimate derivatives \(\tilde{u}'(t)\) from spline
  3. Match ODE: \(\tilde{u}'(t) \approx f(\tilde{u}(t), p, t)\)

Advantages:

  • Very fast (no ODE solving)
  • Good for dense data
  • Useful for initial parameter estimates

Limitations:

  • Less accurate for sparse data
  • Doesn’t account for error accumulation
  • Usually used as first stage before shooting method

When to Use:

  • Dense observational data
  • Quick initial parameter estimates
  • Two-stage method: collocation → shooting

23.8 Key Takeaways

  1. Time-varying confounding requires special methods (G-methods)
  2. G-computation simulates outcomes under treatment strategies
  3. IPTW creates pseudo-populations where treatment is independent of confounders
  4. Challenges (positivity, misspecification, unmeasured confounding) require careful handling
  5. G-methods bridge Observable (data) and Structural (interventions)

23.9 Further Reading

  • Robins (1986): “A new approach to causal inference in mortality studies”
  • Robins et al. (2000): “Marginal structural models”
  • Hernán and Robins (2020): Causal Inference: Comprehensive treatment
  • Chernozhukov et al. (2024): Applied Causal Inference Powered by ML and AI, Chapter 5 (conditional ignorability, propensity scores, and ATE under static confounding; for ML nuisances and DML, see Ch. 9)
  • Rackauckas and Nie (2017): Parameter estimation and inverse problems in differential equations
  • Rackauckas (2026): Parallel Computing and Scientific Machine Learning: comprehensive treatment of inverse problems, multiple shooting, and collocation methods
  • Little and Rubin (2019): Statistical Analysis with Missing Data
  • Laan and Rose (2011): Targeted learning under incomplete observation
  • TMLE and Doubly Robust Estimation: Robust estimation methods
  • Model Validation with Observable Data: Validating inferences
  • Software Patterns: package layout for missingness policies and certificates
  • Incomplete observation: stratum × level treatment of missing, \(R\), and gaps
  • Parameter Estimation via Shooting Method: Shooting method for parameter estimation