26  Model Validation with Observable Data

Status: Draft

v0.4

26.1 Introduction

Model criticism is essential: all models are wrong, but some are useful (Gelman et al. 2013). This chapter develops diagnostics for what is wrong with a fitted model and whether it remains useful for the question at hand. That is part of Seeing in the Observable stratum: validating inferences drawn from actualised data.

26.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?

26.3 Calibration

26.3.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.

26.3.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)

26.3.3 Implementation: Calibration Plots

We can check calibration by comparing predicted probabilities to 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 = 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 26.1: Calibration plot: predicted vs observed probabilities

26.4 Posterior Predictive Checks (PPCs)

26.4.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.

26.4.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

26.4.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

26.4.4 Implementation: Posterior Predictive Checks

We can compare observed data to data simulated from the fitted model:

# 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 = 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)

26.4.5 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.505  (4.294 s)

Figure 26.2: Turing posterior predictive: distribution of replicate means versus the observed mean (vertical line).

Figure 26.3: Posterior predictive check: comparing observed to simulated data

26.5 Residual Analysis

26.5.1 What Are Residuals?

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

26.5.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

26.5.3 Residual Structure

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

Implication: Model is missing important features

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

26.5.4 Implementation: Residual Analysis

We can check residuals for 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 = 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 26.4: Residual analysis: checking for structure in residuals

26.6 Out-of-Domain Validation

26.6.1 The Problem

Models often fail when applied to new domains:

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

26.6.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

26.6.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 26.5: Pooled calibration (mean spectrum OLS, country LOSO).

Figure 26.6: Out-of-domain MAE by held-out country (mean OLS).

Figure 26.7: 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, sec-agescm-tmle), not pooled accuracy alone.

26.7 Stratum context

This chapter addresses Seeing in the Observable stratum: validating what we’ve learned from actualised data. Model validation uses observable data to test whether inferred mechanisms are consistent with observations, bridging the Observable stratum (what we observe) with the Structural/Dynamical strata (what mechanisms exist).

26.8 Key Takeaways

1. Calibration: Predicted probabilities should match observed frequencies
2. Posterior predictive checks: Compare simulated to observed data
3. Residual analysis: Check for structure in residuals
4. Out-of-domain validation: Test generalisation to new domains
5. Model criticism is essential: Find what’s wrong before using the model

26.9 Further Reading

• 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

Gabry, Jonah, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman. 2019. “Visualization in Bayesian Workflow.” Journal of the Royal Statistical Society: Series A 182 (2): 389–402.

Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed. Chapman & Hall/CRC.

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32.