29  Counterfactual Reasoning: Unit-Level Alternatives

Status: Draft

v0.4

29.1 Introduction

Counterfactuals ask: “What would have happened for this specific unit under an alternative?” (Pearl 2009; Imbens and Rubin 2015) Level 3 (Imagining) requires an explicit structural model and a declared cross-world coupling. Shared \(\mathbf u\) is one canonical SCM coupling, not a criterion of identity. Here the Observable role conditions unit-level outcomes on recorded \(Y_t\).

29.2 Counterfactuals as Unit-Level Replays

Technically, counterfactuals compare trajectories for the same unit under different interventions. When we simulate a counterfactual trajectory \(Y^{\iota}_{1:T}(\mathbf{u})\), we replay the same mechanisms with the same \(\mathbf{u}\) but with a different intervention.

29.2.1 Counterfactuals Across Layers

In the three modelling roles (Structural, Dynamical, Measurement), counterfactuals conditioned on records involve:

• For a unit-level query: A declared coupling, often shared \(\mathbf{u}\) (not an identity criterion)
• Under alternative conditions: Different interventions applied to the same mechanisms

The counterfactual trajectory \(Y^{\iota}_{1:T}(\mathbf{u})\) represents the unit’s alternative outcome under intervention \(\iota\), holding \(\mathbf{u}\) fixed.

29.2.2 The Process of Counterfactual Simulation

When we compute a counterfactual:

1. Infer latent state and exogenous history: From observations, infer a posterior over latent states and, when needed, exogenous histories consistent with the data under the model. Do not treat that posterior as a uniquely recovered history.

2. Simulate under an alternative intervention: Under the declared cross-world coupling, use the coupled \(\mathbf{u}\) and a different \(do(\cdot)\) on \(a\) to simulate the alternative trajectory.

3. Compare: Compare the counterfactual trajectory to the observed trajectory (or to another counterfactual), depending on the query.

29.2.3 Why Shared Exogenous Inputs Matter

Fixed \(\mathbf{u}\) is a declared cross-world coupling, not a label for an organism or population. Abduction generally gives a posterior over compatible draws rather than a unique value. The glossary records the limited process interpretation of this construction.

Population-level estimands (TMLE, g-computation) average over many organisms, typically many draws of \(\mathbf{u}\). Counterfactuals fix \(\mathbf{u}\); do not confuse a sampling frame with a shared mechanism class (Chapter 9).

29.3 Counterfactual Trajectories

29.3.1 Definition

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

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

Interpretation: “A declared cross-world coupling, often shared \(\mathbf{u}\), with a different intervention.”

29.3.2 Key Requirement

Counterfactuals require an explicitly declared cross-world coupling. A common SCM coupling holds the exogenous realisation \(\mathbf{u}\) fixed:

• Shared \(\mathbf{u}\) is a modelling convention, not a testable criterion of organismal identity
• Different worlds \(\Leftrightarrow\) different interventions
• Counterfactual: Same \(\mathbf{u}\), different \(do(\cdot)\)

29.4 Computing Counterfactuals

There are two main approaches to computing counterfactuals:

1. Unit-level approach: Infer exogenous inputs \(\mathbf{u}\) for each unit, then simulate a counterfactual for that unit
2. Population-level approach: Estimate population counterfactual means directly from data (e.g., TMLE, g-methods)

29.4.1 Unit-Level Approach

Steps:

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

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

29.4.2 Implementation: Graph Structure and Exogenous Inputs

The causal graph structure determines which exogenous variables must be shared for counterfactual reasoning:

# 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 DAGMakie CairoMakie CausalDynamics Graphs

# Example: Counterfactual requires shared U
# Graph: U → X, U → Y, X → Y
g, _ = confounding_graph(["U", "X", "Y"])  # 1=U, 2=X, 3=Y

# To compute Y^{do(X=0)}(u) given Y^{do(X=1)}(u) = y_obs:
# We need the same u for both worlds

# The graph shows U affects both X and Y
# This means U must be shared across counterfactual worlds

# Check what variables affect Y (these determine which U's must be shared)
parents_Y = get_parents(g, 3)
println("Parents of Y: ", parents_Y)  # Set([1, 2]) = {U, X}

# For counterfactual Y^{do(X=x)}(u), we need:
# - Same U (exogenous noise affecting Y)
# - Different X (intervention changes X)

# The ancestors of Y tell us all variables that could affect Y:
ancestors_Y = get_ancestors(g, 3)
println("Ancestors of Y: ", ancestors_Y)  # Set([1, 2]) = {U, X}

# This shows: for counterfactual reasoning about Y,
# we fix a declared exogenous context U and vary X through intervention

# Package path for additive neural mechanisms on codes:
# MechanismSpec(:generative) → abduce_noise / mechanism_counterfactual
# (CausalDynamics Lux weakdep; see SciML recipes). Raw images encode first.

# Visualise with triangle layout so U → Y is visible
let
    fig, ax, p = dagplot(g;
        figure_size = (600, 400),
        layout = Point2f[Point2f(0.0, 1.0), Point2f(-1.0, 0.0), Point2f(1.0, 0.0)],
        node_color = [RGBf(0.85, 0.65, 0.05), RGBf(0.30, 0.55, 0.70), RGBf(0.30, 0.55, 0.70)],
        labels = ["U", "X", "Y"],
        node_size = 42,
        edge_width = 2.0,
        arrow_size = 16,
        nlabels_fontsize = 18,
        nlabels_color = :white,
        nlabels_align = (:center, :center),
        nlabels_distance = 0,
        padding = 0.40,
    )
    fig  # Only this gets displayed
end

Parents of Y: Set([2, 1])
Ancestors of Y: Set([2, 1])

Figure 29.1: Counterfactual in the Observable stratum: a declared shared-noise coupling supports an L3 comparison

29.4.3 Shared-\(\mathbf{u}\) trajectories on the host-immunity CDM

The same requirement appears on the book’s host-immunity CDM: retain the declared exogenous context and change \(do(A)\). Figure 29.2 shows factual assay \(Y_t\), interventional trajectories under \(do(a=1)\) and \(do(a=0)\), and the shared-\(\mathbf{u}\) counterfactual under treatment (Chapter 28). The shared context specifies the model’s cross-world coupling; it is not a measure of creative advance or an identity criterion.

1:60

Figure 29.2: Host-immunity CDM: factual assay \(Y_t\), interventional trajectories under \(do(a=1)\) and \(do(a=0)\), and a shared-\(\mathbf{u}\) counterfactual under treatment.

29.4.4 Implementation: Graph Structure for Counterfactual Reasoning

The causal graph structure determines what information is needed for counterfactual reasoning. Path and adjustment helpers illustrate that in code:

# 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 CausalDynamics DAGMakie Graphs

# Example: Treatment counterfactual
# Graph: U → A, U → Y, A → Y
g, _ = confounding_graph(["U", "A", "Y"])  # 1=U, 2=A, 3=Y

# To compute counterfactual Y^{do(A=0)}(u) for unit with observed Y^{do(A=1)}(u) = y_obs:
# 1. We need to infer u from observations
# 2. The Markov boundary of Y tells us what variables are needed

mb_Y = markov_boundary(g, 3)  # Outcome Y
println("Markov boundary of Y: ", mb_Y)  # Set([1, 2]) = {U, A}

# This tells us: to reason about Y counterfactually, we need U and A
# If U is unobserved, counterfactuals are not fully identified

# Check if A → Y is identifiable (necessary for counterfactuals)
adj_set = backdoor_adjustment_set(g, 2, 3)
println("Adjustment set for A → Y: ", adj_set)  # Set([1]) = {U}

# If U is unobserved, we cannot identify the causal effect,
# and therefore cannot compute counterfactuals

Markov boundary of Y: Set([2, 1])
Adjustment set for A → Y: Set([1])

29.4.5 Implementation: Available observations and identification

The adjustment query must specify which variables are recorded. In the chunk-counterfactual-identifiability-observable chunk, \(U\) is unobserved and \(L\) is recorded. The set \(\{U,L\}\) blocks the backdoor paths if both variables are available, but adjusting for \(L\) alone leaves \(A \leftarrow U \rightarrow Y\) open. This graph therefore supplies no observed backdoor adjustment set. A unit-level counterfactual requires further structural and cross-world assumptions; an adjustment check alone does not establish its identification (Pearl 2009).

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

g = DiGraph(4)  # 1=U (latent), 2=A, 3=Y, 4=L (recorded)
for (src, dst) in [(1, 2), (1, 3), (2, 3), (4, 2), (4, 3)]
    add_edge!(g, src, dst)
end
observed_covariates = Set([4])
adj_set = backdoor_adjustment_set(g, 2, 3)
# For this graph the only available candidates are {} and {L}.
observed_adjustment_exists = any(
    z -> is_valid_adjustment_set(g, 2, 3, z),
    [Set{Int}(), observed_covariates],
)
@assert adj_set == Set([1, 4])
@assert !observed_adjustment_exists
println("Adjustment set if all covariates were observed: ", adj_set)
println("Backdoor adjustment using recorded covariates: ", observed_adjustment_exists)

Adjustment set if all covariates were observed: Set([1, 4])
Backdoor adjustment using recorded covariates: false

Inferring a posterior over \(U\) requires specified mechanisms, a noise model, and an observation likelihood. It can yield a model-dependent counterfactual distribution. Additional data and assumptions must justify those mechanisms; state-space inference alone does not remove unmeasured confounding.

29.4.6 Implementation: Patient Counterfactual Example

# 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 DAGMakie CairoMakie CausalDynamics Graphs

# Patient counterfactual graph
# Nodes: 1=U (unobserved confounder), 2=A (treatment), 3=Y (outcome), 4=L (observed confounder)
g = DiGraph(4)
add_edge!(g, 1, 2)  # U → A
add_edge!(g, 1, 3)  # U → Y
add_edge!(g, 2, 3)  # A → Y
add_edge!(g, 4, 2)  # L → A
add_edge!(g, 4, 3)  # L → Y

# To compute Y^{do(A=0)}(u) given Y^{do(A=1)}(u) = y_obs:
# 1. Check the recorded adjustment set, L; U remains unavailable.
observed_adjustment_valid = is_valid_adjustment_set(g, 2, 3, Set([4]))
@assert !observed_adjustment_valid
println("Adjustment for recorded L sufficient: ", observed_adjustment_valid)

# 2. Find what variables are needed
adj_set = backdoor_adjustment_set(g, 2, 3)
println("Adjustment set: ", adj_set)  # Set([1, 4]) = {U, L}

# 3. Problem: U is unobserved
# A model-dependent analysis requires an explicit likelihood and prior for U.
# The graph structure shows that U affects both A and Y,
# so we can potentially infer U from (A, Y, L) if we have:
# - Structural assignments f_A(U, L) and f_Y(A, U, L)
# - Assumptions and data that identify or constrain these mechanisms

# 4. Integrate counterfactual outcomes over the posterior of U:
#    - Given: Y^{do(A=1)}(u) = y_obs
#    - Infer: u from y_obs, a_obs=1, l_obs, and known mechanisms
#    - Compute: Y^{do(A=0)}(u) using same u but A=0

# The Markov boundary lists local graph neighbours, not an identification certificate:
mb_Y = markov_boundary(g, 3)
println("Markov boundary of Y: ", mb_Y)  # Set([1, 2, 4]) = {U, A, L}

# Visualise: observed good control L via adjustment colouring; U is latent (unobserved)
let
    fig, ax, p = dagplot(g;
        figure_size = (600, 400),
        layout = Point2f[
            Point2f(-0.8, 1.0),  # U (latent)
            Point2f(-1.0, 0.0),  # A
            Point2f(1.0, 0.0),   # Y
            Point2f(0.8, 1.0),   # L (observed good control)
        ],
        color_by = :adjustment,
        exposure = 2,
        outcome = 3,
        adjustment = Set([4]),
        labels = ["U", "A", "Y", "L"],
        edge_routing = Dict(
            (1, 3) => CurvedEdge(bow = 0.14, side = :right),
            (4, 3) => CurvedEdge(bow = 0.12, side = :left),
        ),
        node_size = 40,
        edge_width = 2.0,
        arrow_size = 15,
        nlabels_fontsize = 16,
        padding = 0.40,
    )
    fig  # Only this gets displayed
end

Adjustment for recorded L sufficient: false
Adjustment set: Set([1, 4])
Markov boundary of Y: Set([2, 1, 4])

Patient counterfactual: observed good control \(L\) for \(A\to Y\); latent \(U\) remains unobserved (Cinelli et al. (Cinelli et al. 2024)).

29.4.7 Implementation: Partial Identification

The graph check in the chunk-partial-identification-observable chunk distinguishes an open causal path from an open backdoor path. It does not calculate numerical bounds; those require an outcome scale and additional assumptions.

# 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 CausalDynamics Graphs

# Example: Counterfactual with unobserved confounder
# Graph: U → A, U → Y, A → Y (U unobserved)
# Nodes: 1=U, 2=A, 3=Y
g = DiGraph(3)
add_edge!(g, 1, 2)  # U → A
add_edge!(g, 1, 3)  # U → Y
add_edge!(g, 2, 3)  # A → Y

# With U unobserved, the empty set is the only available adjustment set.
observed_adjustment_valid = is_valid_adjustment_set(g, 2, 3, Set{Int}())
@assert !observed_adjustment_valid
println("Backdoor adjustment from observed data: ", observed_adjustment_valid)

adj_set = backdoor_adjustment_set(g, 2, 3)
println("Required adjustment: ", adj_set)  # Set([1]) = {U}

# Problem: U is unobserved, so we cannot adjust for it
# Result: Counterfactuals are not fully identified

# However, we can still reason about bounds:
# - The graph structure shows U is a confounder
# - We can use sensitivity analysis to bound counterfactual values
# - The d-separation structure shows what independence assumptions hold

# Check d-separation: A and Y are not d-separated (confounded)
println("A ⫫ Y (no adjustment): ", CausalDynamics.d_separated(g, 2, 3, []))  # false

# Conditioning on U blocks the backdoor, but the direct A → Y remains.
@assert !CausalDynamics.d_separated(g, 2, 3, [1])
println("A ⫫ Y | U in the original graph: ", CausalDynamics.d_separated(g, 2, 3, [1]))
# Remove the outgoing treatment edge only for the backdoor graph check.
g_backdoor = copy(g)
rem_edge!(g_backdoor, 2, 3)
@assert CausalDynamics.d_separated(g_backdoor, 2, 3, [1])

# This structure tells us:
# - Full counterfactuals require U (unobserved)
# - Partial identification is possible via sensitivity analysis
# - Bounds depend on assumptions about U's distribution

Backdoor adjustment from observed data: false
Required adjustment: Set([1])
A ⫫ Y (no adjustment): false
A ⫫ Y | U in the original graph: false

29.5 Partial Identification and Bounds

29.5.1 When Full Identification Isn’t Possible

Sometimes counterfactuals are not fully identified from data. In these cases, we can compute bounds, ranges of possible values.

29.5.2 Counterfactual data as an experimental resource

The usual hierarchy assumes that data are observational or interventional. Some environments permit a stronger design: a counterfactual randomisation can alter what one mechanism receives while leaving another mechanism’s natural value available for conditioning. Such a design produces physically realisable Layer 3 data rather than merely simulating a counterfactual after fitting an SCM.

This distinction yields two useful consequences. First, identifiability can be characterised relative to the counterfactual distributions that the environment can actually generate. Secondly, when a query remains unidentified, those data can shrink its partial-identification bounds. The result is not permission to treat every fitted alternative trajectory as observed evidence: the randomisation must be admissible, its unaffected mechanisms must be stated, and the resulting data regime must be kept distinct from the shared-\(\mathbf{u}\) coupling used for modelled unit-level replay (Raghavan and Bareinboim 2026).

29.5.3 Bounds Methods

• Manski bounds: Worst-case bounds under minimal assumptions
• Sensitivity analysis: How do bounds change with assumptions?
• Instrumental variables: Tighter bounds with additional structure

29.5.4 Interpretation

Bounds tell us what we can know (the range) even when we can’t identify the exact value.

29.6 Stratum context

Given recorded data, structural assumptions, and a declared coupling, counterfactual reasoning determines which unit-level alternative outcomes the model implies.

29.7 Key Takeaways

1. Counterfactuals are declared cross-world replays: coupled exogenous draws, different intervention on \(a\)
2. Not population averages: Distinct from Level-2 estimands that integrate over \(P(\mathbf{u})\)
3. Society vs population: A cohort is a statistical grouping. A Whiteheadian society is more specific: a pattern of continuity carried forward across occasions. Shared membership or a common mechanism alone does not establish that (Ch. 9).
4. Partial identification: Bounds when full identification isn’t possible
5. Counterfactuals bridge data and model: Use data plus assumptions to reason about alternatives

29.8 Further Reading

• Pearl (2009): Causality: Counterfactual framework
• Imbens and Rubin (2015): Causal Inference: Counterfactual methods
• Richardson and Robins (2013): “Single world intervention graphs”
• Hypothesis Generation from Counterfactuals: Using counterfactuals to generate hypotheses
• Experimental Design: Optimal Measurements: Designing studies based on counterfactual reasoning

Cinelli, Carlos, Andrew Forney, and Judea Pearl. 2024. “A Crash Course in Good and Bad Controls.” Sociological Methods & Research 53 (3): 1071–104. https://doi.org/10.1177/00491241221099552.

Imbens, Guido W., and Donald B. Rubin. 2015. Causal Inference in Statistics, Social, and Biomedical Sciences. Cambridge University Press.

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

Raghavan, Arvind, and Elias Bareinboim. 2026. “Causal Identification from Counterfactual Data: Completeness and Bounding Results.” arXiv Preprint arXiv:2602.23541. https://arxiv.org/pdf/2602.23541.

Richardson, Thomas S., and James M. Robins. 2013. “Single World Intervention Graphs (SWIGs): A Unification of the Counterfactual and Graphical Approaches to Causality.” Center for the Statistics and the Social Sciences, University of Washington Series, no. 128.