6  Structural Causal Models as Executable Mechanisms

Status: Draft

v0.6

6.1 Introduction

This chapter develops Structural Causal Models (SCMs) as executable mechanisms: a graph for structure, assignments for how each variable is generated, and exogenous noise for unit-level variation left outside the model (Pearl 2009; Peters et al. 2017). It builds on the graph patterns of Graph Theory and Causal Patterns. In these pages the working type is CausalDynamics.jl’s GraphSCM, with simulate_scm, do_intervention / apply_intervention, and compute_counterfactual for shared-\(U\) unit-level alternatives.

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.

6.2 What Is an SCM?

An SCM is a tuple \(\mathcal{M} = (G, U, F, P(U))\) where (Pearl 2009; Peters et al. 2017):

  • \(G\): Directed acyclic graph (DAG) representing causal structure
  • \(U\): Exogenous (unobserved) variables capturing unmodelled variation and noise
  • \(F\): Structural assignments (functions) encoding the causal mechanisms
  • \(P(U)\): Distribution over exogenous variables

The graph \(G\) encodes which directed dependencies exist; the assignments \(F\) encode how parents and exogenous noise determine each endogenous variable.

6.3 Structural Assignments

Each endogenous variable \(X_i\) is assigned via a structural equation:

\[ X_i \coloneqq f_i(\mathrm{Pa}(X_i), U_i) \]

where \(\mathrm{Pa}(X_i)\) are the parents of \(X_i\) in \(G\). The parent set specifies which variables enter the mechanism; \(f_i\) specifies the mapping.

6.4 Structural Equations as Mechanisms

Each structural equation encodes a causal mechanism from:

  • Parents \(\mathrm{Pa}(X_i)\) (direct causes under the assumed graph)
  • Exogenous noise \(U_i\) (unmodelled variation, stochasticity, or latent inputs)

6.4.1 Defining a graph and mechanisms

The toy DAG \(X \rightarrow Y \leftarrow Z\) has structural equations. In general,

\[ \begin{aligned} X &\coloneqq U_X, \\ Z &\coloneqq U_Z, \\ Y &\coloneqq \beta_{XY}\, X + \beta_{ZY}\, Z + U_Y, \end{aligned} \]

and a numeric specialisation used below is \(\beta_{XY}=2\), \(\beta_{ZY}=3\).

# 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

# Nodes: 1 = X, 2 = Y, 3 = Z
g_fork = SimpleDiGraph(3)
add_edge!(g_fork, 1, 2)  # X → Y
add_edge!(g_fork, 3, 2)  # Z → Y
fork_labels = ["X", "Y", "Z"]

equations_fork = Dict{Int, Function}(
    1 => (u) -> u,                         # X := U_X
    3 => (u) -> u,                         # Z := U_Z
    2 => (x, z, u) -> 2x + 3z + u,         # Y := 2X + 3Z + U_Y (parents sorted)
)
scm_fork = GraphSCM(g_fork, equations_fork, Set{Int}())
U_fork = Dict(1 => 1.0, 2 => 0.1, 3 => 2.0)
vals_fork = simulate_scm(scm_fork, U_fork)
println("Factual values: X=$(vals_fork[1]), Y=$(vals_fork[2]), Z=$(vals_fork[3])")
Factual values: X=1.0, Y=8.1, Z=2.0

Toy SCM \(X \rightarrow Y \leftarrow Z\). Nodes \(X\) and \(Z\) have empty parent sets in \(G\); \(Y\) is determined by both parents and \(U_Y\).

GraphSCM stores a DiGraph, a Dict{Int, Function} of assignments (parent values in index order, then that node’s exogenous draw), and an optional set of exogenous node indices. simulate_scm(scm, U) evaluates in topological order for one fixed realisation of \(U\).

6.5 Modularity Principle

Under intervention \(do(X_i = x)\), the assignment for \(X_i\) is replaced by a constant:

\[ X_i \coloneqq x \quad \text{(instead of } f_i(\mathrm{Pa}(X_i), U_i)\text{)} \]

All other mechanisms remain unchanged. That modularity of causal mechanisms is the local-change semantics of \(do(\cdot)\) (Pearl 2009): the intervened node no longer depends on its natural parents, but descendants still depend on its fixed value.

6.5.1 Implementation: modularity under intervention

NoteFrom CausalDynamics.jl

The modularity principle above is not only prose: apply_intervention is that local change. Incoming edges to the intervened node are removed, and its structural equation becomes a constant assignment. Other equations are left untouched.

# packages/CausalDynamics.jl/src/scm/interventions.jl (excerpt)
function apply_intervention(scm::GraphSCM, intervention::DoIntervention)
    var = intervention.variable
    val = intervention.value
    node = var isa Int ? var : throw(ArgumentError(
        "Symbol-based variable resolution not yet supported for GraphSCM. Use integer node indices.",
    ))

    # Mutilated graph: remove all incoming edges to the intervened node
    new_graph = copy(scm.graph)
    for parent in collect(Graphs.inneighbors(new_graph, node))
        Graphs.rem_edge!(new_graph, parent, node)
    end

    # Replace only this node's equation with a constant
    new_equations = copy(scm.equations)
    new_equations[node] = (args...) -> val

    return GraphSCM(new_graph, new_equations, scm.exogenous)
end

The book’s examples call this API; the excerpt is the mechanism itself. Full reference: CausalDynamics docs.

On the fork, \(do(X = 5)\) removes \(X\)’s natural assignment (and any incoming edges) while \(Y\)’s equation still uses the fixed \(X\) and \(Z\):

scm_do = apply_intervention(scm_fork, do_intervention(1, 5.0))
vals_do = simulate_scm(scm_do, U_fork)
println("Under do(X = 5): X=$(vals_do[1]), Y=$(vals_do[2]), Z=$(vals_do[3])")
println("Parents of X after surgery: ", collect(inneighbors(scm_do.graph, 1)))
println("Y still has parents X and Z: ", has_edge(scm_do.graph, 1, 2) && has_edge(scm_do.graph, 3, 2))
Under do(X = 5): X=5.0, Y=16.1, Z=2.0
Parents of X after surgery: Int64[]
Y still has parents X and Z: true

6.6 Shared-\(U\) counterfactuals

Interventional queries ask what happens in a population under \(do(\cdot)\), typically averaging over \(U\). Unit-level counterfactuals hold the same \(U\) fixed and change only the intervention (Pearl’s abduction–action–prediction; creative advance held constant in the process gloss of Table 8). CausalDynamics exposes that as compute_counterfactual.

The chain \(Z \rightarrow X \rightarrow Y\) makes the arithmetic transparent:

g_chain = DiGraph(3)
add_edge!(g_chain, 1, 2)  # Z → X
add_edge!(g_chain, 2, 3)  # X → Y
equations_chain = Dict{Int, Function}(
    1 => (u) -> u,
    2 => (z, u) -> z + u,
    3 => (x, u) -> 2x + u,
)
scm_chain = GraphSCM(g_chain, equations_chain, Set{Int}())
U_chain = Dict(1 => 1.0, 2 => 0.5, 3 => -0.3)

result = compute_counterfactual(scm_chain, do_intervention(2, 10.0), U_chain)
println("Factual Y = ", result.factual[3])
println("Counterfactual Y under do(X = 10) = ", result.counterfactual[3])
println("Same U_Y in both worlds: exogenous draw for node 3 held at ", U_chain[3])
Factual Y = 2.7
Counterfactual Y under do(X = 10) = 19.7
Same U_Y in both worlds: exogenous draw for node 3 held at -0.3

Factual \(Y = 2(1.0 + 0.5) + (-0.3) = 2.7\); under \(do(X = 10)\) with the same \(U\), \(Y = 2\cdot 10 + (-0.3) = 19.7\). Chapter 7 develops the twin-network reading of this construction.

Chain \(Z \rightarrow X \rightarrow Y\) used for the shared-\(U\) counterfactual. Intervention \(do(X = 10)\) mutilates incoming edges to \(X\) while \(Y\) remains a child of \(X\).

6.7 Linear SCMs and matrix representations

When structural equations are linear, they can be written in matrix form, connecting the graph of Chapter 3 to executable assignments. For a static linear SEM (Peters et al. 2017)

\[ \mathbf{X} = B\mathbf{X} + \mathbf{U}, \]

with \(B_{ij}\) the structural weight on \(X_j \rightarrow X_i\). Under a topological order \(B\) is strictly triangular; its support matches the edge set of \(G\). Time-indexed linear dynamics (Part II) use related sparse transition maps. Matrix analysis and sparse factorisation are the usual algebraic companions (Strang 2016; Horn and Johnson 2013; Golub and Van Loan 2013; Davis 2006).

Sparse dependencies. In most real systems each variable depends on only a small subset of others. Missing edges imply structural zeros in \(B\); present edges mark the only non-zeros that need storing or estimating.

6.7.1 Algebraic path labels, then numbers

The confounding DAG \(X_1 \rightarrow X_2 \rightarrow X_3\), \(X_1 \rightarrow X_3\) is the same bridge as in Chapter 3: first structural path coefficients, then a numeric instance that GraphSCM can simulate.

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 SparseArrays

# Linear SCM: X₁ → X₂ → X₃, X₁ → X₃
g_lin = SimpleDiGraph(3)
add_edge!(g_lin, 1, 2)
add_edge!(g_lin, 2, 3)
add_edge!(g_lin, 1, 3)
labels_lin = ["X₁", "X₂", "X₃"]

# Graphs.edges order: 1→2, 1→3, 2→3
elabels_lin_alg = structural_edge_labels(
    g_lin,
    ["\\beta_{21}", "\\beta_{31}", "\\beta_{32}"];
    latex = true,
)
B_lin_alg = [
    "0"           "0"           "0";
    "\\beta_{21}" "0"           "0";
    "\\beta_{31}" "\\beta_{32}" "0";
]

B_lin = [
    0.0  0.0  0.0;   # X₁ := U₁
    0.5  0.0  0.0;   # X₂ := 0.5 X₁ + U₂
    0.3  0.7  0.0;   # X₃ := 0.3 X₁ + 0.7 X₂ + U₃
]

equations_lin = Dict{Int, Function}(
    1 => (u) -> u,
    2 => (x1, u) -> 0.5 * x1 + u,
    3 => (x1, x2, u) -> 0.3 * x1 + 0.7 * x2 + u,
)
scm_lin = GraphSCM(g_lin, equations_lin, Set{Int}())
U_lin = Dict(1 => 1.0, 2 => 0.0, 3 => 0.0)
vals_lin = simulate_scm(scm_lin, U_lin)

A = adjacency_matrix(g_lin)
println("Graphs.jl adjacency A (A[i,j]=1 ⇒ i→j):")
println(Matrix(A))
println("\nCoefficient matrix B:")
println(B_lin)
println("Non-zeros in sparse(B): ", nnz(sparse(B_lin)), " of ", length(B_lin))
println("simulate_scm with U=(1,0,0): X₁=$(vals_lin[1]), X₂=$(vals_lin[2]), X₃=$(vals_lin[3])")
println("  (check: X₂ = 0.5, X₃ = 0.3 + 0.7*0.5 = 0.65)")

"""Unshaded matrix grid; bold non-zeros / non-\"0\" entries (SEM B defaults)."""
function coefficient_matrix_grid!(
    ax, B, tick_labels;
    digits = 1,
    fontsize = 18,
    xlabel = "parent j",
    ylabel = "child i",
)
    n = size(B, 1)
    for i in 1:n, j in 1:n
        poly!(ax, Rect(j - 0.5, i - 0.5, 1, 1);
            color = :white, strokecolor = :black, strokewidth = 1.25)
        entry = B[i, j]
        if entry isa Real
            text_str = string(round(entry; digits = digits))
            is_zero = iszero(entry)
            label = text_str
        else
            text_str = String(entry)
            is_zero = text_str == "0" || text_str == "0.0"
            label = occursin(r"[\\^_{}]", text_str) ? Makie.latexstring(text_str) : text_str
        end
        text!(ax, j, i;
            text = label,
            align = (:center, :center),
            color = :black,
            fontsize = fontsize,
            font = is_zero ? :regular : :bold,
        )
    end
    ax.xticks = (1:n, tick_labels)
    ax.yticks = (1:n, tick_labels)
    ax.xlabel = xlabel
    ax.ylabel = ylabel
    ax.yreversed = true
    ax.aspect = DataAspect()
    xlims!(ax, 0.5, n + 0.5)
    ylims!(ax, 0.5, n + 0.5)
    return ax
end
Graphs.jl adjacency A (A[i,j]=1 ⇒ i→j):
[0 1 1; 0 0 1; 0 0 0]

Coefficient matrix B:
[0.0 0.0 0.0; 0.5 0.0 0.0; 0.3 0.7 0.0]
Non-zeros in sparse(B): 3 of 9
simulate_scm with U=(1,0,0): X₁=1.0, X₂=0.5, X₃=0.6499999999999999
  (check: X₂ = 0.5, X₃ = 0.3 + 0.7*0.5 = 0.65)
Main.Notebook.coefficient_matrix_grid!

General linear confounding / mediation triangle: structural path coefficients on the DAG, symbolic \(B\), and assignments. \(B_{ij}\neq 0\) means \(X_j\rightarrow X_i\) (Peters et al. 2017).

Numeric instance \(\beta_{21}=0.5\), \(\beta_{31}=0.3\), \(\beta_{32}=0.7\), executed by GraphSCM / simulate_scm for \(U=(1,0,0)\).

6.7.2 Executable SCMs at scale: GraphPPL and latent interfaces

For high-dimensional observations (e.g. MIRS spectra Y), the SCM is still defined over structural nodes, but the inference layer should not duplicate thousands of wavenumbers inside a PPL. One workable pattern (AgeSCM; Chapters 28b and 30):

  1. Identification / graph: DAG and adjustment sets (CausalDynamics helpers such as prepare_for_rxinfer, prepare_for_tmle).
  2. Encoder: e.g. Flux ẑ = g(Y; nuisances residualised).
  3. Probabilistic programme: GraphPPL / RxInfer @model on ẑ and the outcome of interest, not on raw Y (live wiring in Chapter 30; poultry ATE in Chapter 5).

Interventions remain structural; the encoder is a learned observation mechanism whose outputs feed the PPL.

6.8 Stratum context

This chapter addresses doing in the Structural layer: interventions as local changes to a data-generating process while other mechanisms stay fixed. Identification of what can be learned from data is the subject of Chapter 5; time-indexed CDMs extend the same modular idea across occasions in Part II and Chapter 28.

6.9 Summary

An SCM pairs a DAG with structural assignments and exogenous noise. GraphSCM makes that tuple executable: simulate_scm settles one unit for fixed \(U\), apply_intervention implements modularity by mutilating parents of the intervened node, and compute_counterfactual compares factual and alternative outcomes under shared \(U\). Linear maps \(\mathbf{X}=B\mathbf{X}+\mathbf{U}\) mirror the adjacency pattern with structural path coefficients on the edges; at high dimension, keep the PPL on a low-dimensional structural interface rather than raw high-dimensional measurements.

6.10 Further Reading