21  Counterfactual Dynamics: Alternative Trajectories

Status: Draft

v0.4

21.1 Introduction

This chapter explores counterfactual reasoning in dynamical systems: alternative concrescences for the same organism (fixed creative advance \(\mathbf{u}\)) under different interventions on \(A_t\) (Chapter 9). This is Imagining in the Dynamical stratum: alternative trajectories of \(\mu_{1:T}\), not population averages.

21.2 Counterfactual Trajectories in Dynamical Systems

21.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: Same organism (same creative advance \(\mathbf{u}\)), different intervention on \(a\), same society \((f, h, G)\).

21.2.2 Key Requirement

Counterfactuals in dynamical systems require shared creative advance \(\mathbf{u}\):

  • Same organism \(\Leftrightarrow\) same \(\mathbf{u}\)
  • Different intervention \(\Leftrightarrow\) different \(do(\cdot)\) on \(a\)
  • Same society \(\Leftrightarrow\) same structural/dynamical mechanisms

21.3 Computing Counterfactual Dynamics

21.3.1 Unit-Level Approach

Steps:

  1. Infer \(\mathbf{u}\): From observed trajectory, infer the exogenous noise realisation for this unit
  2. Simulate counterfactual: With same \(\mathbf{u}\), simulate under alternative intervention
  3. Compare trajectories: Compare counterfactual trajectory to observed trajectory

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

21.3.2 Implementation: Computing Counterfactual Trajectories

Here’s an example showing how to compute 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 21.1: Counterfactual dynamics: comparing observed and counterfactual trajectories for a specific unit

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

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

21.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 21.2: Bounds for counterfactual dynamics when full identification isn’t possible

21.5 Stratum context

This chapter addresses Imagining in the Dynamical stratum: what alternative dynamic trajectories are possible? Counterfactual dynamics explore what would have happened for specific units under different interventions, given the same exogenous noise realisation.

21.6 Key Takeaways

  1. Counterfactual trajectories: Alternative dynamic trajectories for specific units
  2. Creative advance \(\mathbf{u}\): Same organism across alternative concrescences
  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

21.7 Further Reading