25  Model Validation with Observable Data

Status: Draft

v0.4

25.1 Introduction

All models are wrong, but some remain adequate for a stated question (Gelman et al. 2013). This chapter develops diagnostics for what a fitted model misses and for whether those failures matter to the intended inference. That is Seeing in the Observable stratum: validating inferences drawn from records.

The same discipline applies whether the fitted object is a TMLE/LMTP nuisance, a state-space filter, a Turing or RxInfer head, or a learned CDM. Applied practice is iterative: start with a small generative story, check it against data and against data you simulated yourself, expand when checks fail, and re-check (Gelman et al. 2026; Gabry et al. 2019). The sections below are tools for that loop, not a one-pass checklist after a single fit.

25.2 Why Model Criticism?

Models are approximations. Model criticism helps us:

  1. Find failures: Where does the model break?
  2. Assess usefulness: Is the model good enough for the question?
  3. Guide improvement: What should we fix?
  4. Build trust: Can we trust the model’s predictions?

25.2.1 Criticism as iteration in a CDM

A default order for CDMs and their observable heads is:

  1. Simulate from the current structure (known mechanisms, or a candidate \(f\) and \(h\)) and recover parameters or estimands you already know (oracle \(\tau\), known \(\sigma_w\), known shift effect).
  2. Push forward predictions or posterior draws before conditioning on the real \(Y\) (prior or prior-predictive style checks when the head is Bayesian).
  3. Fit on recorded data, then run calibration, residual, and posterior predictive checks (below).
  4. Expand the graph, the dynamics, or the observation map only when a check fails for a scientifically decisive feature.
  5. Re-check, including out-of-domain splits when transport is part of the claim.

That loop is the observable half of building a CDM. Identification (Part I) still decides which queries are even estimable; criticism decides whether the fitted mechanisms are good enough for those queries (Gelman et al. 2026). The loop is nested rather than sealed: a failed check can send the analysis back to a graph, estimand, temporal declaration, intervention, or observation map.

25.2.2 What a failed check tells us

Different checks support different claims. A successful check is evidence about the component it targets, not a general certificate for the causal model as a whole.

Failure or instability Immediate interpretation What it does not establish
Poor calibration or posterior predictive fit The observation model, residual structure, or missing measurement process may be inadequate That the causal graph is wrong by itself
Poor oracle or parameter recovery The simulator, estimator, or implementation may not reproduce the declared mechanism That the field data are uninformative
Failed simulation-based calibration The inference implementation, posterior computation, or parameterisation may be at fault That the scientific mechanism is false
Poor held-out trajectory or perturbation response The transition law, latent-state representation, or initial-state model may be inadequate That the estimand is unidentified
Weak positivity or unstable policy curve The treatment support, policy range, estimand, or design may not sustain the intended contrast That predictive accuracy is poor
Unstable result under confounding or graph alternatives The identification assumptions, graph, or sensitivity range may be doing the work That a single fitted model is computationally invalid
Poor out-of-domain performance The transport assumptions or domain-varying mechanisms may not hold That the source-domain estimand is unidentified

The distinction matters because model criticism, computational validation, causal identification, and transport assessment answer different questions, even when the evidence from one bears on the others.

25.3 Simulated-data checks and calibration of generators

Before trusting criticism of real assays, ask whether the pipeline can recover truth when truth is known. Simulate a CDM or a simple observational DGP with a known interventional contrast, run the same identification and estimation path you intend for the field data, and compare the estimate to the oracle. Failures here are software or model-class bugs, not mysteries in the cohort.

When the head is Bayesian (Turing, RxInfer, or a hierarchical observation model), simulation-based calibration goes further: draw parameters from the prior, simulate data, refit, and check that posterior quantiles are uniformly calibrated under repeated draws (Gelman et al. 2026). You need not run a full SBC battery in every chapter example; the habit matters. Stress notebooks under CausalTargeted already use oracle recovery for continuous MTP; the same idea applies to PPL heads and to UDE training on trajectories generated from a known CDM (Ch. 28).

These checks sit beside positivity and sensitivity reports for LMTP (Ch. 19, 23). They answer “does this generative + inferential stack behave?” rather than “are backdoors blocked?”

25.3.1 Oracle recovery of a known backdoor effect

The confounded Gaussian DGP \(Z\to X\), \(Z\to Y\), \(X\to Y\) with structural slope \(\tau=2\) is the same golden path as in Chapter 30. Simulate observational rows, run CausalDynamics infer_backdoor_effect, and compare \(\hat\tau\) to the known mechanism. Large absolute error is a pipeline failure.

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

Random.seed!(21)
n = 500
τ_oracle = 2.0
z = randn(n)
x = z .+ 0.4 .* randn(n)
y = τ_oracle .* x .+ z .+ 0.3 .* randn(n)
df_oracle = DataFrame(Z = z, X = x, Y = y)
g_oracle = DiGraph(3)
add_edge!(g_oracle, 1, 2)
add_edge!(g_oracle, 1, 3)
add_edge!(g_oracle, 2, 3)
nm_oracle = Dict(1 => :Z, 2 => :X, 3 => :Y)

t_oracle = @elapsed begin
    global rec = infer_backdoor_effect(
        g_oracle, df_oracle, 2, 3;
        node_names = nm_oracle, iterations = 30,
    )
end
abs_err = abs(rec.τ_mean - τ_oracle)
println(
    "Oracle recovery: identifiable = ", rec.identifiable,
    "  confounders = ", rec.confounders,
    "  τ̂ = ", round(rec.τ_mean; digits = 3),
    "  |τ̂ − τ| = ", round(abs_err; digits = 3),
    "  (", round(t_oracle; digits = 3), " s)",
)
Oracle recovery: identifiable = true  confounders = [:Z]  τ̂ = 1.998  |τ̂ − τ| = 0.002  (8.339 s)
Figure 25.1: Oracle recovery: variational backdoor mean \(\hat\tau\) against the structural slope \(\tau=2\) used to simulate the table.

25.4 Calibration

25.4.1 What Is Calibration?

A model is calibrated if predicted probabilities match observed frequencies (Gelman et al. 2013).

Example: If the model predicts 80% probability of an event, the event should occur ~80% of the time.

25.4.2 How to Check

  • Calibration plots: Predicted vs observed probabilities (Gelman et al. 2013)
  • Reliability diagrams: Binned predictions vs observed frequencies
  • Scoring rules: Brier score, log score (Vehtari et al. 2017)

25.4.3 Implementation: Calibration Plots

Calibration is assessed by comparing predicted probabilities with observed frequencies:

# 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(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))

@auto_using Random Distributions GLM DataFrames CairoMakie

Random.seed!(42)

# Simulate data
n = 1000
X = rand(Normal(0, 1), n)
p_true = 1 ./ (1 .+ exp.(-(0.5 .+ 0.8 .* X)))  # True probability
Y = rand.(Bernoulli.(p_true))

df = DataFrame(X = X, Y = Y)

# Fit model
model = glm(@formula(Y ~ X), df, Binomial(), LogitLink())
p_pred = GLM.predict(model, df)  # Predicted probabilities

# Bin predictions and compute observed frequencies
n_bins = 10
bins = collect(range(0, 1, length = n_bins + 1))
bin_centers = (bins[1:end-1] .+ bins[2:end]) ./ 2
observed_freq = Float64[]
predicted_mean = Float64[]

for i in 1:n_bins
    # Handle last bin to include upper bound
    if i == n_bins
        in_bin = (p_pred .>= bins[i]) .& (p_pred .<= bins[i+1])
    else
        in_bin = (p_pred .>= bins[i]) .& (p_pred .< bins[i+1])
    end
    if sum(in_bin) > 0
        push!(observed_freq, mean(Y[in_bin]))
        push!(predicted_mean, mean(p_pred[in_bin]))
    end
end

# Compute Brier score
brier_score = mean((p_pred .- Y).^2)
println("Calibration assessment:")
println("  Brier score: ", round(brier_score, digits=4))
println("  (Lower is better; perfect = 0)")
Calibration assessment:
  Brier score: 0.2089
  (Lower is better; perfect = 0)
Figure 25.2: Calibration plot: predicted vs observed probabilities

25.5 Posterior Predictive Checks (PPCs)

25.5.1 The Idea

Compare observed data to data simulated from the fitted model (Gelman et al. 2013; Gabry et al. 2019): \[ Y^{\text{rep}} \sim P(Y \mid \theta_{\text{posterior}}) \]

If \(Y^{\text{rep}}\) looks like \(Y^{\text{obs}}\), the model captures the data.

25.5.2 Test Statistics

Choose test statistics \(T(Y)\) that capture important features:

  • Means, variances: First and second moments
  • Autocorrelations: Temporal structure
  • Extrema: Rare events
  • Domain-specific: Scientific quantities of interest

25.5.3 Interpretation

  • \(p\)-values: Extreme values (near 0 or 1) indicate problems
  • Visual comparison: Do simulated data look like observed data?
  • Multiple statistics: Check many aspects

25.5.4 Implementation: Posterior Predictive Checks

The fitted model is checked by comparing observed data with data simulated from it:

# 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 observed data
n = 200
X = rand(Normal(0, 1), n)
Y_true = 1.0 .+ 0.5 .* X .+ rand(Normal(0, 0.3), n)

df = DataFrame(X = X, Y = Y_true)

# Fit model
model = lm(@formula(Y ~ X), df)
Y_pred = GLM.predict(model, df)
σ_residual = std(residuals(model))

# Simulate from fitted model (posterior predictive)
n_sims = 1000
Y_rep = zeros(n, n_sims)

for sim in 1:n_sims
    Y_rep[:, sim] = Y_pred .+ rand(Normal(0, σ_residual), n)
end

# Test statistics
function test_stat_mean(Y)
    """Test statistic: mean of observations."""
    return mean(Y)
end

function test_stat_max(Y)
    """Test statistic: maximum of observations."""
    return maximum(Y)
end

T_obs_mean = test_stat_mean(Y_true)
T_obs_max = test_stat_max(Y_true)

T_rep_mean = [test_stat_mean(Y_rep[:, sim]) for sim in 1:n_sims]
T_rep_max = [test_stat_max(Y_rep[:, sim]) for sim in 1:n_sims]

# p-values (proportion of T_rep >= T_obs)
p_value_mean = mean(T_rep_mean .>= T_obs_mean)
p_value_max = mean(T_rep_max .>= T_obs_max)

println("Posterior predictive check:")
println("  Mean: p-value = ", round(p_value_mean, digits=3))
println("  Maximum: p-value = ", round(p_value_max, digits=3))
println("  (Extreme p-values near 0 or 1 indicate problems)")
Posterior predictive check:
  Mean: p-value = 0.499
  Maximum: p-value = 0.707
  (Extreme p-values near 0 or 1 indicate problems)

25.5.5 Expand when a check fails

The iterate-and-expand habit needs a failing check, not only a green one. Here \(Y = 1 + 0.8\,X + \varepsilon\) is simulated once. An intercept-only working model generates replicates that cannot reproduce the observed association with \(X\); the test statistic \(T = |\mathrm{cor}(Y,X)|\) is extreme. Expanding to \(Y\sim X\) restores agreement. That is the same loop as in the Introduction of this chapter: criticise, expand, re-check (Gelman et al. 2026).

# Reuse the linear DGP pattern; independent seed for a clean expand demo
using Statistics
Random.seed!(7)
n_exp = 200
X_exp = rand(Normal(0, 1), n_exp)
Y_exp = 1.0 .+ 0.8 .* X_exp .+ rand(Normal(0, 0.35), n_exp)
df_exp = DataFrame(X = X_exp, Y = Y_exp)
T_obs_cor = abs(cor(Y_exp, X_exp))
n_sims_exp = 800

function ppc_cor_p(model, df; n_sims = n_sims_exp)
    # PPC p-value for |cor(Y, X)| under Gaussian replicates from `model`.
    # Use `GLM.predict(model)` (original frame). `GLM.predict(model, df)` fails for
    # intercept-only formulae: StatsModels builds an empty term NamedTuple.
    Ŷ = GLM.predict(model)
    σ̂ = std(residuals(model))
    T_rep = Float64[]
    for _ in 1:n_sims
        Y_rep = Ŷ .+ rand(Normal(0, σ̂), nrow(df))
        push!(T_rep, abs(cor(Y_rep, df.X)))
    end
    return mean(T_rep .>= T_obs_cor), T_rep
end

model_wrong = lm(@formula(Y ~ 1), df_exp)
model_ok = lm(@formula(Y ~ X), df_exp)
p_wrong, T_rep_wrong = ppc_cor_p(model_wrong, df_exp)
p_ok, T_rep_ok = ppc_cor_p(model_ok, df_exp)

println("PPC for |cor(Y, X)| (observed T = ", round(T_obs_cor; digits = 3), "):")
println("  Intercept-only: p = ", round(p_wrong; digits = 3), "  (expect near 0 → expand)")
println("  Y ~ X:          p = ", round(p_ok; digits = 3), "  (expect mid-range → retain)")
PPC for |cor(Y, X)| (observed T = 0.924):
  Intercept-only: p = 0.0  (expect near 0 → expand)
  Y ~ X:          p = 0.532  (expect mid-range → retain)
Figure 25.3: Expand when a check fails: replicate distribution of \(|\mathrm{cor}(Y,X)|\) under an intercept-only model (left) versus \(Y\sim X\) (right). Amber line: observed association.

25.5.6 Turing PPC (same linear story)

A short NUTS fit yields posterior draws of \((\mu_0,\beta,\sigma)\) for the same DGP; replicate datasets are then draws from the posterior predictive:

using Turing: Turing, @model, sample, NUTS
using Distributions
using Statistics

Turing.@model function lin_reg(X, Y)
    μ0 ~ Normal(0, 2)
    β ~ Normal(0, 2)
    σ ~ truncated(Normal(0.5, 0.5); lower = 0.05)
    for i in eachindex(Y)
        Y[i] ~ Normal(μ0 + β * X[i], σ)
    end
end

t_ppc = @elapsed begin
    global ch_ppc = sample(lin_reg(X, Y_true), NUTS(0.8), 300; progress = false)
end
μ0s = vec(Array(ch_ppc[:μ0]))
βs = vec(Array(ch_ppc[:β]))
σs = vec(Array(ch_ppc[:σ]))
n_rep = 200
T_rep_mean_t = Float64[]
for k in 1:n_rep
    i = rand(1:length(μ0s))
    Yrep = μ0s[i] .+ βs[i] .* X .+ σs[i] .* randn(n)
    push!(T_rep_mean_t, mean(Yrep))
end
p_t = mean(T_rep_mean_t .>= mean(Y_true))
println("Turing PPC mean p-value = ", round(p_t; digits = 3), "  (", round(t_ppc; digits = 3), " s)")
Turing PPC mean p-value = 0.55  (4.246 s)
Figure 25.4: Turing posterior predictive: distribution of replicate means versus the observed mean (vertical line).
Figure 25.5: Posterior predictive check: comparing observed to simulated data

25.6 Residual Analysis

25.6.1 What Are Residuals?

Residuals measure how well the model fits: \[ r_t = Y_t - \mathbb{E}[Y_t \mid X_t, \theta] \]

25.6.2 What to Check

  • Mean: Should be near zero
  • Variance: Should be constant (homoscedasticity)
  • Autocorrelation: Should be near zero (no temporal structure)
  • Distribution: Should match observation model

25.6.3 Residual Structure

Problem: Residuals show structure (trends, cycles, correlations)

Implication: Model is missing important features

Solution: Add missing components (trends, seasonality, interactions)

25.6.4 Implementation: Residual Analysis

Residuals can reveal structure that indicates model problems:

# 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 missing quadratic term (model misspecification)
n = 200
X = rand(Uniform(-2, 2), n)
Y = 1.0 .+ 0.5 .* X .+ 0.3 .* X.^2 .+ rand(Normal(0, 0.2), n)  # True model has X²

df = DataFrame(X = X, Y = Y)

# Fit misspecified model (missing X² term)
model = lm(@formula(Y ~ X), df)
Y_pred = GLM.predict(model, df)
residuals = Y .- Y_pred

# Check residuals
residual_mean = mean(residuals)
residual_std = std(residuals)

# Check autocorrelation (for time series, here we check against X)
# In time series, check lag-1 autocorrelation
# Here, we check correlation with X to detect missing terms
residual_X_corr = cor(residuals, X)

println("Residual analysis:")
println("  Mean: ", round(residual_mean, digits=4), " (should be ≈ 0)")
println("  Std: ", round(residual_std, digits=4))
println("  Correlation with X: ", round(residual_X_corr, digits=4), " (should be ≈ 0)")
println("\nInterpretation:")
if abs(residual_X_corr) > 0.1
    println("  ⚠️  Residuals correlated with X → missing term (e.g., X²)")
else
    println("  ✓ Residuals appear independent of X")
end
Residual analysis:
  Mean: -0.0 (should be ≈ 0)
  Std: 0.4089
  Correlation with X: 0.0 (should be ≈ 0)

Interpretation:
  ✓ Residuals appear independent of X
Figure 25.6: Residual analysis: checking for structure in residuals

25.7 Out-of-Domain Validation

25.7.1 The Problem

Models often fail when applied to new domains:

  • Different time periods: Temporal generalisation
  • Different populations: Population generalisation
  • Different conditions: Condition generalisation

Predictive generalisation is not causal identification. A model can forecast well in a new country and still answer the wrong interventional query if adjustment or transport assumptions fail; conversely, an identified estimand can be estimated poorly under shift. Out-of-domain checks (below) probe the former; graphs, certificates, and positivity/sensitivity reports probe the latter (Gelman et al. 2026; Pearl 2009). AgeSCM pairs both: LOSO error and IPTW overlap.

25.7.2 How to Validate

  • Temporal split: Train on past, validate on future
  • Population split: Train on one population, validate on another
  • Condition split: Train on one condition, validate on another

25.7.3 AgeSCM: calibration and out-of-domain checks

Case Study 3 validates age models on held-out countries (T, S, B, M), plus replica and experiment splits. Phase 4 of the AgeSCM repository reports:

  • Calibration: reliability of predicted vs observed age (mean OLS and Stage 1), pooled and per domain.
  • OOD MAE: bar charts of hold-out domain error for country LOSO.
  • Residuals: pooled residual structure for mean OLS on country split.
Figure 25.7: Pooled calibration (mean spectrum OLS, country LOSO).
Figure 25.8: Out-of-domain MAE by held-out country (mean OLS).
Figure 25.9: Residuals (mean OLS, country LOSO).

Interpretation: Models are not well calibrated to ±2 days across all domains; S is easiest, M (small n) volatile. Transport claims must pair these diagnostics with IPTW overlap tables (Ch 20, Section 24.3.4), not pooled accuracy alone.

25.8 Stratum context

Model validation uses observable data to test whether inferred mechanisms are consistent with observations. Predictive fit can expose a model failure, but it does not establish causal structure on its own.

25.9 Key Takeaways

  1. Iterate: Simulate and recover, fit, criticise, expand, re-check (see oracle recovery and the intercept-only versus \(Y\sim X\) PPC)
  2. Calibration: Predicted probabilities should match observed frequencies
  3. Posterior predictive checks: Compare simulated to observed data
  4. Residual analysis: Check for structure in residuals
  5. Out-of-domain validation: Test predictive generalisation; keep it distinct from identification
  6. Model criticism is essential: Find what is wrong before using the model for \(do(\cdot)\) or policy claims
  7. Interpret checks narrowly: Each check bears on the component and claim it was designed to test; none is a blanket validation of the graph

25.10 Further Reading

  • Gelman et al. (2026): Bayesian Workflow: iterative building, checking, simulated-data experimentation (practice applies beyond Bayesian heads)
  • Gelman et al. (2013): Bayesian Data Analysis: Model criticism
  • Gabry et al. (2019): “Visualization in Bayesian workflow”
  • Vehtari et al. (2017): “Practical Bayesian model evaluation”
  • Observational Methods: Learning from Data: Methods that need validation
  • TMLE and Doubly Robust Estimation: Robust methods that still need validation