27  Interventional Reasoning: Forecasting Under Interventions

Status: Draft

v0.4

27.1 Learning Objectives

After reading this chapter, you will be able to:

• Distinguish conditional forecasting from interventional forecasting
• Forecast under interventions in dynamical systems
• Understand how structural changes propagate through time
• Recognise when conditional and interventional forecasts differ

27.2 Introduction

This chapter focuses on Doing in the Observable world—how can we forecast what will happen under interventions? A critical distinction in causal-dynamical systems: conditional forecasting (what happens if I observe X?) vs interventional forecasting (what happens if I do X?) (Robins 1986; Robins et al. 2000; Pearl 2009). This chapter demonstrates the operational difference and provides templates for interventional queries.

27.3 Conditional vs Interventional Forecasting

27.3.1 Conditional Forecasting

Conditional forecasting asks: “What will happen if I observe X?”

\[ P(Y_{t+1} \mid Y_{1:t}, A_t = a) \]

Interpretation: Distribution of future outcomes given past observations and observed treatment value.

When to use: - Prediction: Forecast future observations - No intervention: We’re not changing anything, just observing - Association: Level 1 of causal hierarchy

27.3.2 Interventional Forecasting

Interventional forecasting asks: “What will happen if I do X?”

\[ P^{do(A_t = a)}(Y_{t+1} \mid Y_{1:t}) \]

Interpretation: Distribution of future outcomes under intervention \(do(A_t = a)\).

When to use: - Policy evaluation: What happens if we implement a policy? - Treatment effects: What happens if we give treatment? - Intervention: Level 2 of causal hierarchy

27.3.3 When Do They Differ?

If treatment assignment is confounded, conditional and interventional forecasts differ (Pearl 2009; Robins et al. 2000):

• Conditional: \(P(Y \mid A = a)\) adjusts for confounders via conditioning
• Interventional: \(P^{do(A = a)}(Y)\) sets treatment structurally, breaking confounder links

Key insight: Conditional forecasting observes what happens naturally; interventional forecasting predicts what happens when we change the system.

27.3.4 Implementation: Conditional vs Interventional Forecasting

Here’s an example demonstrating the difference when treatment is confounded:

# 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 confounding
# L: Disease severity (confounder)
# A: Treatment (depends on L - confounding!)
# Y: Outcome (depends on L and A)

n = 1000
L = rand(Normal(0, 1), n)  # Confounder
# Treatment assignment: more severe (higher L) → more likely to treat
p_treat = 1 ./ (1 .+ exp.(-(L .- 0.5)))
A = rand.(Bernoulli.(p_treat))
# Outcome: treatment helps, but severity also affects outcome
Y = 0.5 .* A .- 0.8 .* L .+ rand(Normal(0, 0.2), n)

df = DataFrame(L = L, A = A, Y = Y)

# Conditional forecast: P(Y | A = 1)
# This adjusts for confounders via conditioning
conditional_model = lm(@formula(Y ~ A + L), df)
Y_cond_A1 = predict(conditional_model, DataFrame(A = ones(n), L = L))
conditional_forecast = mean(Y_cond_A1)

# Interventional forecast: P^{do(A = 1)}(Y)
# This sets A = 1 structurally, breaking the L → A link
# We need to marginalise over L: E[Y | do(A=1)] = E_L[E[Y | A=1, L]]
# Use marginal distribution of L (not conditional on A=1)
# Sample L from its marginal distribution
L_marginal = rand(Normal(mean(L), std(L)), n)  # Marginal distribution of L
Y_interv_A1 = predict(conditional_model, DataFrame(A = ones(n), L = L_marginal))
interventional_forecast = mean(Y_interv_A1)

# The forecasts differ because:
# - Conditional: P(Y | A=1) conditions on observed A=1, which is correlated with L
# - Interventional: P^{do(A=1)}(Y) sets A=1, breaking the correlation

# Visualise the difference
let
    fig = Figure(size = (800, 400))
    ax = Axis(fig[1, 1], title = "Conditional vs Interventional Forecasts", 
              xlabel = "Method", ylabel = "Expected Outcome")
    
    barplot!(ax, [1, 2], [conditional_forecast, interventional_forecast], 
             label = ["Conditional\nP(Y|A=1)", "Interventional\nP^{do(A=1)}(Y)"], 
             color = [:blue, :red])
    axislegend(ax)
    
    fig  # Only this gets displayed
end

println("Forecast comparison:")
println("  Conditional: E[Y | A=1] = ", round(conditional_forecast, digits=3))
println("  Interventional: E[Y | do(A=1)] = ", round(interventional_forecast, digits=3))
println("\nNote: When treatment is confounded, these differ!")
println("  Conditional adjusts for confounders via conditioning")
println("  Interventional breaks the confounder link structurally")

Forecast comparison:
  Conditional: E[Y | A=1] = 0.539
  Interventional: E[Y | do(A=1)] = 0.547

Note: When treatment is confounded, these differ!
  Conditional adjusts for confounders via conditioning
  Interventional breaks the confounder link structurally

27.4 Interventional Forecasting in Dynamical Systems

27.4.1 Structural Changes Propagate Through Time

When we intervene \(do(A_t = a)\), we modify the structural equation for \(A_t\):

\[ A_t \coloneqq a \quad \text{(instead of } \pi(H_t, U^A_t) \text{)} \]

This structural change propagates forward through time:

1. Immediate effect: \(A_t\) is set to \(a\)
2. Direct effect: \(A_t\) affects \(Y_t\) and future states
3. Indirect effects: Changes propagate through feedback loops
4. Long-term effects: System evolves toward new attractor

27.4.2 Computing Interventional Forecasts

Steps: 1. Set intervention: \(do(A_t = a)\) for all \(t\) in forecast period 2. Simulate forward: Use state transition \(f\) to propagate effects 3. Account for feedback: Past outcomes affect future states 4. Aggregate: Compute expected outcomes or distributions

27.4.3 Example: Treatment Intervention

Forecast patient outcome under treatment intervention:

\[ P^{do(\text{Treatment}_t = \text{drug A})}(\text{Recovery}_{t+1} \mid \text{History}_{1:t}) \]

This requires: - Understanding how treatment affects recovery (causal mechanism) - Accounting for time-varying confounding - Propagating effects through time

27.4.4 Implementation: Interventional Forecasting in Dynamical Systems

We can forecast under interventions by simulating the system forward with the intervention applied:

# 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 DifferentialEquations CairoMakie

# Example: Disease model with treatment intervention
# Forecast: What happens if we intervene with treatment at t=5?

β = 0.2  # Disease progression
α = 0.3  # Treatment effectiveness

# Original system (no intervention)
function disease_original!(du, u, p, t)
    """Disease model: original dynamics with confounding (treatment depends on severity)."""
    X = u[1]
    # Treatment assignment depends on severity (confounding)
    A = X > 0.5 ? 1.0 : 0.0
    du[1] = β * X - α * A * X
end

# System under intervention do(A = 1) for t >= 5
function disease_intervened!(du, u, p, t)
    """Disease model: intervened dynamics with do(A = 1) starting at t=5."""
    X = u[1]
    # Intervention: do(A = 1) starting at t=5
    A = t >= 5.0 ? 1.0 : (X > 0.5 ? 1.0 : 0.0)
    du[1] = β * X - α * A * X
end

u0 = [0.3]  # Initial severity
tspan = (0.0, 20.0)

# Solve both systems
prob_original = ODEProblem(disease_original!, u0, tspan)
sol_original = solve(prob_original, Tsit5())

prob_intervened = ODEProblem(disease_intervened!, u0, tspan)
sol_intervened = solve(prob_intervened, Tsit5())

# Extract final-time forecasts for reuse (for both plotting and summary text)
final_original = sol_original.u[end][1]
final_intervened = sol_intervened.u[end][1]

# Visualise forecasts
let
    fig = Figure(size = (1000, 400))
    ax1 = Axis(fig[1, 1], title = "Disease Severity Over Time", 
               xlabel = "Time", ylabel = "Severity")
    ax2 = Axis(fig[1, 2], title = "Forecast Comparison at t=20", 
               xlabel = "Method", ylabel = "Severity")
    
    lines!(ax1, sol_original.t, [u[1] for u in sol_original.u], 
           label = "Original (conditional)", linewidth = 2, color = :blue)
    lines!(ax1, sol_intervened.t, [u[1] for u in sol_intervened.u], 
           label = "Intervened (do(A=1) at t=5)", linewidth = 2, color = :red, linestyle = :dash)
    vlines!(ax1, [5.0], color = :gray, linestyle = :dot, linewidth = 1, label = "Intervention")
    axislegend(ax1)
    
    barplot!(ax2, [1, 2], [final_original, final_intervened], 
             label = ["Conditional", "Interventional"], color = [:blue, :red])
    axislegend(ax2)
    
    fig  # Only this gets displayed
end

println("Interventional forecast:")
println("  Severity at t=20 (original): ", round(final_original, digits=3))
println("  Severity at t=20 (intervened): ", round(final_intervened, digits=3))
println("  Improvement: ", round((final_original - final_intervened) / final_original * 100, digits=1), "%")

Interventional forecast:
  Severity at t=20 (original): 0.502
  Severity at t=20 (intervened): 0.106
  Improvement: 78.9%

27.5 Uncertainty in Interventional Forecasts

When making interventional forecasts, we must account for uncertainty (Rackauckas 2026). This uncertainty comes from:

1. Parameter uncertainty: Uncertainty about model parameters
2. Model uncertainty: Uncertainty about model structure
3. Measurement uncertainty: Error in observations used to fit the model
4. Numerical uncertainty: ODE solver errors

Why it matters: Intervention effect estimates without uncertainty bounds are incomplete. We need to know not just the expected effect, but also how uncertain we are about it.

27.5.1 Prediction Intervals for Interventions

Prediction intervals quantify uncertainty in intervention effects:

\[ P^{do(A=a)}(Y_{t+1} \in [L, U] \mid Y_{1:t}) = 1 - \alpha \]

where \([L, U]\) is the \((1-\alpha)\) prediction interval.

Methods: - Monte Carlo: Sample parameters from posterior, run model, compute quantiles - Linearisation: Approximate using first-order Taylor expansion - Bootstrap: Resample data, refit model, compute intervals

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 DifferentialEquations Distributions Random Statistics 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

# Parameter uncertainty: β ~ N(0.3, 0.05²), γ ~ N(0.1, 0.02²)
Random.seed!(42)
n_samples = 1000
β_samples = rand(Normal(0.3, 0.05), n_samples)
γ_samples = rand(Normal(0.1, 0.02), n_samples)

u0 = [0.99, 0.01, 0.0]
tspan = (0.0, 50.0)

# Intervention: reduce transmission rate (vaccination)
β_intervened = 0.15  # 50% reduction

# Monte Carlo: sample parameters, run model, collect outcomes
# Store pairs to ensure matching samples
outcome_pairs = Tuple{Float64, Float64}[]

for i in 1:n_samples
    # Original (no intervention)
    p_orig = (β_samples[i], γ_samples[i])
    prob_orig = ODEProblem(sir!, u0, tspan, p_orig)
    sol_orig = solve(prob_orig, Tsit5())
    
    # Intervened (vaccination)
    p_int = (β_intervened, γ_samples[i])  # β reduced, γ unchanged
    prob_int = ODEProblem(sir!, u0, tspan, p_int)
    sol_int = solve(prob_int, Tsit5())
    
    # Only keep pairs where both succeeded
    if SciMLBase.successful_retcode(sol_orig.retcode) && SciMLBase.successful_retcode(sol_int.retcode) && 
       length(sol_orig.u) > 0 && length(sol_int.u) > 0
        peak_I_orig = maximum([u[2] for u in sol_orig.u])
        peak_I_int = maximum([u[2] for u in sol_int.u])
        if isfinite(peak_I_orig) && isfinite(peak_I_int)
            push!(outcome_pairs, (peak_I_orig, peak_I_int))
        end
    end
end

# Extract separate arrays
outcomes_original = [p[1] for p in outcome_pairs]
outcomes_intervened = [p[2] for p in outcome_pairs]

# Check if we have enough samples
if length(outcomes_original) == 0
    error("No successful ODE solves. Check model parameters and solver settings.")
end

# Compute prediction intervals (95%)
quantiles_orig = quantile(outcomes_original, [0.025, 0.5, 0.975])
quantiles_int = quantile(outcomes_intervened, [0.025, 0.5, 0.975])

# Intervention effect with uncertainty (now guaranteed to have matching pairs)
effect_samples = outcomes_original .- outcomes_intervened
effect_quantiles = quantile(effect_samples, [0.025, 0.5, 0.975])

# Visualise
let
fig = Figure(size = (1200, 400))
ax1 = Axis(fig[1, 1], title = "Peak Infection: Original vs Intervened", 
           xlabel = "Scenario", ylabel = "Peak Infection Rate",
           xticks = ([1, 2], ["Original", "Intervened"]))
ax2 = Axis(fig[1, 2], title = "Intervention Effect Distribution", 
           xlabel = "Reduction in Peak Infection", ylabel = "Density")
ax3 = Axis(fig[1, 3], title = "Prediction Intervals", 
           xlabel = "Scenario", ylabel = "Peak Infection Rate",
           xticks = ([1, 2], ["Original", "Intervened"]))

# Box plots
boxplot!(ax1, [1], outcomes_original, color = :blue, width = 0.3)
boxplot!(ax1, [2], outcomes_intervened, color = :red, width = 0.3)

# Effect distribution
hist!(ax2, effect_samples, bins = 50, normalization = :pdf, color = (:green, 0.5))
vlines!(ax2, [effect_quantiles[2]], color = :black, linestyle = :dash, linewidth = 2, label = "Median")
vlines!(ax2, effect_quantiles[[1, 3]], color = :gray, linestyle = :dot, linewidth = 1, label = "95% CI")

# Prediction intervals
errorbars!(ax3, [1], [quantiles_orig[2]], 
           [quantiles_orig[2] - quantiles_orig[1]], [quantiles_orig[3] - quantiles_orig[2]],
           color = :blue, linewidth = 2)
scatter!(ax3, [1], [quantiles_orig[2]], color = :blue, markersize = 10)

errorbars!(ax3, [2], [quantiles_int[2]], 
           [quantiles_int[2] - quantiles_int[1]], [quantiles_int[3] - quantiles_int[2]],
           color = :red, linewidth = 2)
scatter!(ax3, [2], [quantiles_int[2]], color = :red, markersize = 10)

axislegend(ax2)

fig
end

println("Uncertainty in interventional forecasts:")
println("  Original peak I: ", round(quantiles_orig[2], digits=3), 
        " (95% CI: [", round(quantiles_orig[1], digits=3), ", ", round(quantiles_orig[3], digits=3), "])")
println("  Intervened peak I: ", round(quantiles_int[2], digits=3), 
        " (95% CI: [", round(quantiles_int[1], digits=3), ", ", round(quantiles_int[3], digits=3), "])")
println("  Intervention effect: ", round(effect_quantiles[2], digits=3), 
        " (95% CI: [", round(effect_quantiles[1], digits=3), ", ", round(effect_quantiles[3], digits=3), "])")
println("  Relative reduction: ", round(effect_quantiles[2] / quantiles_orig[2] * 100, digits=1), "%")

Uncertainty in interventional forecasts:
  Original peak I: 0.298 (95% CI: [0.121, 0.496])
  Intervened peak I: 0.061 (95% CI: [0.011, 0.222])
  Intervention effect: 0.226 (95% CI: [0.087, 0.328])
  Relative reduction: 75.9%

27.5.2 Decision-Making Under Uncertainty

When making decisions based on interventional forecasts, we need to account for uncertainty:

• Robust decisions: Choose interventions that work well across parameter uncertainty
• Risk assessment: Consider worst-case scenarios
• Value of information: How much would better parameter estimates improve decisions?

Key Points: - Uncertainty quantification is essential for reliable interventional forecasts - Prediction intervals provide decision-relevant information - Parameter uncertainty propagates through models to create forecast uncertainty - Connects to counterfactual reasoning (see Counterfactual Reasoning: Unit-Level Alternatives)

27.6 World Context

This chapter addresses Doing in the Observable world—how can we forecast under interventions? Interventional forecasting requires understanding how inner worlds (Structural/Dynamical) become outer worlds (Observable), and how we can modify that process through interventions.

27.7 Key Takeaways

1. Conditional vs interventional: Observing vs doing—fundamentally different
2. Interventional forecasting: Predict what happens under interventions
3. Structural changes propagate: Interventions affect system evolution through time
4. Confounding matters: Conditional and interventional forecasts differ when confounded

27.8 Further Reading

• Pearl (2009): Causality — Do-calculus and interventions
• Robins (1986): “A new approach to causal inference”
• Robins et al. (2000): “Marginal structural models”
• Rackauckas (2026): Parallel Computing and Scientific Machine Learning — comprehensive treatment of uncertainty quantification and prediction intervals
• Uncertainty Propagation: Methods for propagating uncertainty through models
• Policy Evaluation and Dynamic Treatment Strategies: Evaluating treatment strategies
• Causal Decision-Making: Using forecasts for decisions

Pearl, Judea. 2009. Causality: Models, Reasoning, and Inference. 2nd ed. Cambridge University Press.

Rackauckas, Christopher. 2026. Parallel Computing and Scientific Machine Learning (SciML): Methods and Applications. https://doi.org/10.5281/zenodo.6917234.

Robins, James M. 1986. “A New Approach to Causal Inference in Mortality Studies with a Sustained Exposure Period—Application to Control of the Healthy Worker Survivor Effect.” Mathematical Modelling 7 (9–12): 1393–512.

Robins, James M., Miguel A. Hernán, and Babette Brumback. 2000. “Marginal Structural Models and Causal Inference in Epidemiology.” Epidemiology 11 (5): 550–60.