20  Counterfactual Dynamics: Alternative Trajectories

Status: Draft

v0.4

20.1 Introduction

Counterfactual trajectories require an explicitly declared cross-world coupling. Shared exogenous draws are the canonical SCM coupling used in the examples, but they are a modelling assumption rather than an empirical criterion of organismal identity. The analysis concerns coupled alternative trajectories, not population averages.

20.2 Counterfactual Trajectories in Dynamical Systems

20.2.1 Definition

For a fixed exogenous realisation \(\mathbf{u}\), the counterfactual trajectory under intervention \(\iota\) in a dynamical system is:

\[ Y^{\iota}_{1:T}(\mathbf{u}) \]

Interpretation: the factual and counterfactual trajectories share a declared coupling, undergo different interventions, and are evaluated under the same declared model version \((G,F,h)\).

20.2.2 Key Requirement

Counterfactuals require abduction, action and prediction:

\[ p(Y^\iota_{t+1:T}\mid y_{0:t}, a_{0:t}, c) = \int p_{\mathcal M^\iota}(Y^\iota_{t+1:T}\mid u)\, p_{\mathcal M}(u\mid y_{0:t},a_{0:t},c)\,du. \]

Shared \(\mathbf u\) is one cross-world coupling. It neither establishes personal identity nor follows uniquely from observed data; alternative admissible couplings should be reported as sensitivity analyses.

20.3 Computing Counterfactual Dynamics

20.3.1 Unit-Level Approach

Steps:

  1. Abduction: Infer a posterior distribution over latent state and exogenous history from the observed trajectory; do not assume that a posterior is a uniquely recovered history
  2. Action: Transform \(\mathcal M\) into \(\mathcal M^\iota\) by the declared intervention
  3. Prediction: Simulate or integrate counterfactual trajectories under the declared coupling
  4. Compare trajectories: Compare posterior counterfactual and factual summaries

Challenge: Inferring \(\mathbf{u}\) from observations requires strong assumptions (full structural model with known noise structure).

20.3.2 Implementation: Computing Counterfactual Trajectories

The following example computes counterfactual trajectories for a specific unit:

# 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 OrdinaryDiffEq Random CairoMakie

Random.seed!(123)

# Example: Treatment timing counterfactual
# Observed: Treatment started at t=10
# Counterfactual: What if treatment started at t=5?

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

# Step 1: Simulate observed trajectory (treatment at t=10)
function disease_observed!(du, u, p, t)
    """Disease model with observed treatment trajectory (treatment at t=10)."""
    X = u[1]
    # Treatment starts at t=10
    A = t >= 10.0 ? 1.0 : 0.0
    du[1] = β * X - α * A * X
end

u0_observed = [0.1]  # Initial severity
tspan = (0.0, 20.0)
prob_observed = ODEProblem(disease_observed!, u0_observed, tspan)
sol_observed = solve(prob_observed, Tsit5())

# Step 2: Infer the "unit" (in this deterministic case, unit = initial condition)
# For deterministic systems, the unit is fully determined by initial conditions
# For stochastic systems, we'd need to infer the noise realisation

# Step 3: Simulate counterfactual (same unit, different intervention: treatment at t=5)
function disease_counterfactual!(du, u, p, t)
    """Disease model with counterfactual treatment trajectory (treatment at t=5)."""
    X = u[1]
    # Counterfactual: Treatment starts at t=5 (earlier)
    A = t >= 5.0 ? 1.0 : 0.0
    du[1] = β * X - α * A * X
end

# Same initial condition (same unit)
u0_counterfactual = u0_observed
prob_counterfactual = ODEProblem(disease_counterfactual!, u0_counterfactual, tspan)
sol_counterfactual = solve(prob_counterfactual, Tsit5())

# Compare outcomes
final_observed = sol_observed.u[end][1]
final_counterfactual = sol_counterfactual.u[end][1]
println("Final severity:")
println("  Observed (treatment at t=10): ", round(final_observed, digits=3))
println("  Counterfactual (treatment at t=5): ", round(final_counterfactual, digits=3))
println("  Improvement from earlier treatment: ", round((final_observed - final_counterfactual) / final_observed * 100, digits=1), "%")
Final severity:
  Observed (treatment at t=10): 0.264
  Counterfactual (treatment at t=5): 0.064
  Improvement from earlier treatment: 75.9%
Figure 20.1: Counterfactual dynamics: comparing observed and counterfactual trajectories for a specific unit

20.3.3 Example: “What Would Have Happened If We Had Intervened Earlier?”

Question: For a specific patient, what would have happened if we had started treatment at time \(t_0\) instead of \(t_1\)?

Process:

  1. Infer \(\mathbf{u}\) from observed trajectory (treatment started at \(t_1\))
  2. Simulate counterfactual: same \(\mathbf{u}\), but treatment starts at \(t_0\)
  3. Compare: counterfactual trajectory vs observed trajectory

20.4 Bounds and Partial Identification

When full counterfactual dynamics are not identified, we can still obtain bounds:

  • Non-parametric bounds: Range of possible counterfactual trajectories
  • Sensitivity parameters: How results change with assumptions about unobserved confounders
  • Uncertainty propagation: How uncertainty in \(\mathbf{u}\) affects counterfactual trajectories

20.4.1 Implementation: Bounds for Counterfactual Dynamics

When we cannot fully identify counterfactual trajectories, we can compute bounds:

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

# Example: Counterfactual with uncertainty in treatment effect
# We know treatment helps, but exact effect α is uncertain: α ∈ [0.2, 0.4]

β = 0.2  # Disease progression (known)

# Lower bound: α = 0.2 (weakest treatment effect)
function disease_lower!(du, u, p, t)
    """Disease model: lower bound (α = 0.2) for partial identification."""
    X = u[1]
    A = t >= 5.0 ? 1.0 : 0.0
    du[1] = β * X - 0.2 * A * X  # Lower bound on treatment effect
end

# Upper bound: α = 0.4 (strongest treatment effect)
function disease_upper!(du, u, p, t)
    """Disease model: upper bound (α = 0.4) for partial identification."""
    X = u[1]
    A = t >= 5.0 ? 1.0 : 0.0
    du[1] = β * X - 0.4 * A * X  # Upper bound on treatment effect
end

u0 = [0.1]
tspan = (0.0, 20.0)

prob_lower = ODEProblem(disease_lower!, u0, tspan)
sol_lower = solve(prob_lower, Tsit5())

prob_upper = ODEProblem(disease_upper!, u0, tspan)
sol_upper = solve(prob_upper, Tsit5())

println("Counterfactual bounds (final severity):")
println("  Lower bound (weakest treatment): ", round(sol_lower.u[end][1], digits=3))
println("  Upper bound (strongest treatment): ", round(sol_upper.u[end][1], digits=3))
println("  Range: ", round(sol_upper.u[end][1] - sol_lower.u[end][1], digits=3))
Counterfactual bounds (final severity):
  Lower bound (weakest treatment): 0.26
  Upper bound (strongest treatment): 0.013
  Range: -0.247
Figure 20.2: Bounds for counterfactual dynamics when full identification isn’t possible

20.5 Stratum context

Counterfactual dynamics describe alternative trajectories for specific units under different interventions and a declared coupling of latent and exogenous histories. Shared values are one convention, not an identity criterion.

20.6 Key Takeaways

  1. Counterfactual trajectories: Alternative dynamic trajectories for specific units
  2. Cross-world coupling: declared shared-noise assumptions support alternative outcomes
  3. Unit-level reasoning: Counterfactuals are for specific units, not populations
  4. Bounds: Partial identification when full identification isn’t possible
  5. Counterfactual dynamics bridge Dynamical and Observable: Use dynamic models to reason about alternative observable outcomes

20.7 Further Reading