10  Counterfactuals: Unit-Level Alternatives at Structural Level

Status: Draft

v0.6

10.1 Introduction

Counterfactuals ask what would have happened for this unit under an alternative intervention (Pearl 2009; Imbens and Rubin 2015). That is Pearl’s third rung: not “what happens on average if we do \(a\)?” (rung 2), but “given what we saw for unit \(i\), what would \(Y\) have been under \(do(a')\)?”

At the Structural level the machinery is an SCM \((G,F,P(\mathbf{U}))\) plus Pearl’s three-step procedure: abduction (recover the unit’s exogenous realisation \(\mathbf{u}\)), action (apply \(do(\cdot)\)), prediction (re-simulate with the same \(\mathbf{u}\)). CausalDynamics.jl packages the last two steps in compute_counterfactual once \(\mathbf{u}\) is in hand; Chapter 4 already showed a short chain example. This chapter slows that process down, with graphs and arithmetic at each step.

NotePackages implement the examples

This book’s framing is general: structure, dynamics, and observation for complex dynamical systems. Executable examples use a Julia stack, CausalDynamics.jl (identify), CausalTargeted.jl (estimate LMTP / mediation), DAGMakie.jl (display). The book is not package documentation; manuals and changelogs live in the repositories.

Code keeps standard Pearl / targeted-learning names (Intervention, do_surgery, backdoor_adjustment_set, run_lmtp_grid, …). Prose may use a sparse process gloss where it helps: prehension (taking account of), shared \(U\) as creative advance (Introduction), L3 alternative concrescences, and Chapter 9’s organism / environment / Markov blanket / society vocabulary for CDM scope (\(\mu\), \(\eta\), \(s\), \(a\), \(G\), \(F\)). Prefer ordinary words (observation, intervention, time step) elsewhere. Mapping: Concept Reference, Tables 4 and 8.

Package source as notation. Where the implementation is the explanation (graph surgery, path blocking, shared-\(U\) counterfactuals), chapters may show a short, syntax-highlighted excerpt from the owned packages in a “From …jl” callout, alongside executed API chunks. Those excerpts are curated under snippets/package-source/; they are not a substitute for the Documenter manuals.

10.2 Interventions versus counterfactuals

Intervention (rung 2) averages over units:

\[ \mathbb{E}[Y \mid do(A=a)] = \int Y^{do(A=a)}(\mathbf{u})\, P(\mathbf{u})\, d\mathbf{u}. \]

Counterfactual (rung 3) holds one unit’s \(\mathbf{u}\) fixed (Richardson and Robins 2013):

\[ Y^{do(A=a)}(\mathbf{u}) \quad\text{versus}\quad Y^{do(A=a')}(\mathbf{u}). \]

The same structural map \(F\) and graph \(G\) appear in both queries; only whether \(\mathbf{u}\) is integrated out or held fixed changes the estimand. Process gloss: fixed \(\mathbf{u}\) is the creative advance of one organism; different \(do(\cdot)\) yields an alternative concrescence (Table 8).

10.3 Pearl’s three steps

Pearl (Pearl 2009, Ch.~7) factors every structural counterfactual into:

  1. Abduction — Infer \(\mathbf{u}\) (or a posterior over \(\mathbf{u}\)) from evidence in the factual world.
  2. Action — Replace the assignment for the intervened variable(s) by a constant; mutilate incoming edges (modularity).
  3. Prediction — Evaluate the mutilated SCM with the same \(\mathbf{u}\).
NoteFrom CausalDynamics.jl

Structural counterfactuals are three lines once abduction has fixed \(\mathbf{u}\): simulate the factual world, mutilate the graph under \(do(\cdot)\), then simulate again with the same exogenous values.

# packages/CausalDynamics.jl/src/scm/counterfactuals.jl (excerpt)
function compute_counterfactual(scm::GraphSCM, intervention::DoIntervention,
                                 exogenous_values::Dict{Int, T}) where T
    factual = simulate_scm(scm, exogenous_values)
    scm_cf = counterfactual_graph(scm, intervention)
    counterfactual = simulate_scm(scm_cf, exogenous_values)
    return (factual = factual, counterfactual = counterfactual)
end

Shared exogenous_values is the creative advance held fixed; only the intervened mechanism changes. Discrete-time CDMs use the same idea via counterfactual(cdm, factual.noise; …) (Chapter 28).

The rest of the chapter walks a single ecological unit through these three steps, then notes what fails when confounders stay latent.

10.4 Worked example: one site, one “what if”

Use the rainfall–herbivore–richness DAG from the interventions practical (Modern Inference): rainfall \(R\) affects herbivore pressure \(H\) and species richness \(S\), and \(H\) also affects \(S\):

\[ \begin{aligned} R &\coloneqq U_R, \\ H &\coloneqq a R + U_H, \\ S &\coloneqq b R + c H + U_S, \end{aligned} \]

with numeric weights \(a=0.5\), \(b=0.7\), \(c=0.4\). Site A is observed at \(R=0.5\), \(H=1\), \(S=1.5\). The counterfactual question is: what would Site A’s richness have been if herbivore pressure had been doubled (\(do(H=2)\))?

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using DAGMakie CairoMakie CausalDynamics Graphs

# Nodes: 1=R, 2=H, 3=S
g = SimpleDiGraph(3)
add_edge!(g, 1, 2)  # R → H
add_edge!(g, 1, 3)  # R → S
add_edge!(g, 2, 3)  # H → S
labels = ["R", "H", "S"]
layout = Point2f[Point2f(0.0, 1.0), Point2f(-1.0, 0.0), Point2f(1.0, 0.0)]

a, b, c = 0.5, 0.7, 0.4
equations = Dict{Int, Function}(
    1 => (u) -> u,
    2 => (r, u) -> a * r + u,
    3 => (r, h, u) -> b * r + c * h + u,
)
scm = GraphSCM(g, equations, Set{Int}())

# Observed Site A (factual world)
R_obs, H_obs, S_obs = 0.5, 1.0, 1.5
(0.5, 1.0, 1.5)

Factual DAG for Site A: rainfall \(R\), herbivores \(H\), richness \(S\). Structural path \(H\rightarrow S\) is the channel we will intervene on.

10.4.1 Step 1: Abduction

Solve the structural equations for the exogenous draws consistent with the observation. Because \(R\) has no parents, \(U_R = R_{\mathrm{obs}}\). Then

\[ \begin{aligned} U_H &= H_{\mathrm{obs}} - a R_{\mathrm{obs}} = 1 - 0.5\cdot 0.5 = 0.75, \\ U_S &= S_{\mathrm{obs}} - b R_{\mathrm{obs}} - c H_{\mathrm{obs}} = 1.5 - 0.7\cdot 0.5 - 0.4\cdot 1 = 0.75. \end{aligned} \]

U_R = R_obs
U_H = H_obs - a * R_obs
U_S = S_obs - b * R_obs - c * H_obs
U = Dict(1 => U_R, 2 => U_H, 3 => U_S)

# Sanity: replaying the factual SCM with this U recovers the observation
factual_check = simulate_scm(scm, U)
@assert factual_check[1] ≈ R_obs
@assert factual_check[2] ≈ H_obs
@assert factual_check[3] ≈ S_obs
println("Abduced U = (U_R=$(U_R), U_H=$(U_H), U_S=$(U_S))")
println("Factual replay: R=$(factual_check[1]), H=$(factual_check[2]), S=$(factual_check[3])")
Abduced U = (U_R=0.5, U_H=0.75, U_S=0.7499999999999999)
Factual replay: R=0.5, H=1.0, S=1.5

Abduction recovers unit-specific exogenous draws. \(U_R\), \(U_H\), and \(U_S\) are fixed for Site A and will be reused under \(do(H=2)\).

For a fully deterministic invertible SCM, abduction is this algebraic solve. With stochastic mechanisms one samples (or approximates) \(P(\mathbf{u}\mid \text{evidence})\). The same Site A observation under a soft Gaussian likelihood yields a posterior concentrated near \((U_H,U_S)=(0.75,0.75)\):

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

Turing.@model function abduct_site(R_obs, H_obs, S_obs, a, b, c)
    U_H ~ Normal(0, 1)
    U_S ~ Normal(0, 1)
    σ = 0.05
    H_obs ~ Normal(a * R_obs + U_H, σ)
    S_obs ~ Normal(b * R_obs + c * H_obs + U_S, σ)
end

t_abd = @elapsed begin
    global ch_abd = sample(abduct_site(R_obs, H_obs, S_obs, a, b, c), NUTS(0.8), 400; progress = false)
end
U_H_post = mean(ch_abd[:U_H])
U_S_post = mean(ch_abd[:U_S])
println("Posterior means: U_H≈", round(U_H_post; digits = 3),
    ", U_S≈", round(U_S_post; digits = 3), "  (", round(t_abd; digits = 3), " s)")
println("Algebraic: U_H=", U_H, ", U_S=", U_S)
Posterior means: U_H≈0.749, U_S≈0.75  (3.735 s)
Algebraic: U_H=0.75, U_S=0.7499999999999999

Stochastic abduction: joint posterior draws of \((U_H,U_S)\) for Site A under a soft Gaussian likelihood. The algebraic solution (cross) sits at the mode.

Push posterior draws of \(\mathbf{u}\) through action and prediction to get a counterfactual distribution, not a single number. For Quarto speed we keep \(400\) NUTS draws (\(\approx 10\,\mathrm{ms}\) after warmup in local benchmarks).

10.4.2 Step 2: Action

Replace the assignment for \(H\) by the constant \(2\) and delete edges into \(H\) (here \(R\rightarrow H\)). Rainfall and the abduced noises are unchanged; only the herbivore mechanism is modularly replaced (Pearl 2009).

ι = do_intervention(2, 2.0)  # do(H = 2)
scm_cf = counterfactual_graph(scm, ι)
println("Parents of H after surgery: ", collect(inneighbors(scm_cf.graph, 2)))
println("S still has parents R and H: ", has_edge(scm_cf.graph, 1, 3) && has_edge(scm_cf.graph, 2, 3))
Parents of H after surgery: Int64[]
S still has parents R and H: true

Action: factual DAG (left) versus mutilated DAG under \(do(H=2)\) (right). Incoming arrows to \(H\) are cut; \(S\) still depends on \(R\) and the fixed \(H\).

10.4.3 Step 3: Prediction

Simulate the mutilated model with the same \(\mathbf{u}\):

\[ S_{\mathrm{cf}} = b R + c\cdot 2 + U_S = 0.7\cdot 0.5 + 0.4\cdot 2 + 0.75 = 1.9. \]

Equivalently, \(S_{\mathrm{cf}} = S_{\mathrm{obs}} + c(2 - H_{\mathrm{obs}}) = 1.5 + 0.4 = 1.9\): doubling herbivores adds one unit of \(c\) to richness for this site.

result = compute_counterfactual(scm, ι, U)
S_fact = result.factual[3]
S_cf = result.counterfactual[3]
println("Factual S = ", S_fact)
println("Counterfactual S under do(H=2) = ", S_cf)
println("Individual effect Δ = ", S_cf - S_fact)
@assert S_cf ≈ 1.9
Factual S = 1.5
Counterfactual S under do(H=2) = 1.9
Individual effect Δ = 0.3999999999999999

Prediction for Site A: observed richness versus counterfactual richness under \(do(H=2)\), same abduced \(\mathbf{u}\). The gap is the unit-level effect \(c\,(2-H_{\mathrm{obs}})=0.4\).

A population regression coefficient for \(H\) reports an average effect. The counterfactual above is tailored to Site A’s abduced \((U_H,U_S)\); another site with different residuals would move by the same structural \(c\cdot\Delta H\) only if mechanisms are linear and identical, but its levels would differ.

10.5 Shared \(\mathbf{u}\) is the unit

Nothing in the action or prediction steps redraws noise. If one sampled a fresh \(\mathbf{u}'\sim P(\mathbf{U})\) after intervening, the query would be interventional (rung 2), not counterfactual (rung 3). In code, that discipline is literally the shared exogenous_values dictionary passed to both simulate_scm calls inside compute_counterfactual.

10.6 When the graph blocks full identification

Abduction needs enough structure to recover (or constrain) \(\mathbf{u}\). If a latent confounder \(U\) sits on a backdoor between treatment and outcome and is never observed, unit-level counterfactuals are typically only partially identified: one can still report bounds or sensitivity parameters rather than a point value (Pearl 2009; Imbens and Rubin 2015).

# Latent U confounds A and Y
g_lat = SimpleDiGraph(3)
add_edge!(g_lat, 1, 2)  # U → A
add_edge!(g_lat, 1, 3)  # U → Y
add_edge!(g_lat, 2, 3)  # A → Y

adj = backdoor_adjustment_set(g_lat, 2, 3)
println("Minimal backdoor set for A→Y: ", adj)  # {U}
println("A ⫫ Y with empty conditioning? ", CausalDynamics.d_separated(g_lat, 2, 3, Int[]))
# Given U, the backdoor is blocked, but the causal path A→Y remains open — so A ⫫̸ Y | U.
# The identification problem is that U is latent: we cannot form that adjustment set from data.
println("U observed in the graph? (node 1 is the confounder we treat as latent in prose)")
Minimal backdoor set for A→Y: Set([1])
A ⫫ Y with empty conditioning? false
U observed in the graph? (node 1 is the confounder we treat as latent in prose)

Latent confounder \(U\) on the backdoor \(A\leftarrow U\rightarrow Y\). The minimal backdoor set is \(\{U\}\); when \(U\) is unobserved, neither population adjustment nor unit-level abduction of \(U\) is available from \((A,Y)\) alone.

Practical posture when \(U\) is missing: state the identifying assumptions, report bounds or E-value-style sensitivity where appropriate, and avoid pretending a point counterfactual is identified from the observational margin alone.

10.7 Stratum context

This chapter is Imagining in the Structural layer: unit-level alternatives implied by an SCM. Related material elsewhere:

  • Dynamical counterfactuals (Chapter 16): alternative trajectories under shared process noise
  • Observable counterfactuals (Chapter 25): prediction targets under estimated mechanisms
  • Twin-network / SWIG diagrams appear again in identification and design chapters when potential outcomes must be drawn explicitly (Richardson and Robins 2013)

10.8 Summary

Structural counterfactuals are interventional queries conditioned on a recovered exogenous realisation. Decompose them as abduction (infer \(\mathbf{u}\)), action (modular \(do(\cdot)\) surgery), and prediction (shared-\(\mathbf{u}\) simulation). compute_counterfactual automates action and prediction once abduction is done; graphs that leave essential noise latent force bounds or sensitivity instead of point answers.

10.9 Further Reading