37  Software Engineering Patterns for Causal-Dynamical Modelling

Status: Draft

v0.4

37.1 Introduction

This appendix covers layout and practice: how the book’s Julia stack is organised, how a thin application wires identification to estimation, and how tests and dependency pins keep work reproducible. It is not a substitute for the conceptual chapters, and it is not the package manuals. Principles transfer to other languages; the concrete layout uses Julia.

37.2 Working stack

In these pages the working stack is CausalDynamics.jl (graphs, identification, SCMs / CDMs), CausalTargeted.jl (cross-fitted LMTP and interventional mediation), and DAGMakie.jl (DAG figures). Exported names match the literature (identify, do_intervention, run_lmtp_grid, dagplot, …). Process vocabulary belongs in this book (Chapter 9; Table 8); package manuals stay Pearl / targeted-learning facing (CausalTargeted methods).

Division of labour:

Layer Owns Does not own
Identification / CDM simulation (CausalDynamics) Graphs, identify, temporal unrolling, GraphSCM / DiscreteTimeCDM Cross-fitted nuisance grids
Estimation (CausalTargeted) LMTP / mediation grids, positivity, MC stability, sensitivity helpers Cohort loaders, DAG string parsers
Application repos Data, registries, concordance Re-implementing EIF estimators

When extending the book or the packages, keep API names stable, add process gloss only in book prose where it clarifies, and update Table 8 if you introduce new central mappings. For manuals and changelogs, use the repositories, not this chapter as API documentation.

TipSmall-n checklist (CausalTargeted)

For continuous MTP and mediation with tens to low hundreds of units (DΓ­az et al. 2023; Laan and Rose 2011):

  1. Start from recommend_run_options(n; engine, n_mediators) (lean Super Learner when \(n < 80\); parallel=false by default).
  2. Inspect positivity / support (positivity_report or grid positivity=true) (Petersen et al. 2012).
  3. For mediation, sweep nested MC (mediation_n_mc_sweep) until signs and SEs stabilise (Liu et al. 2024).
  4. Report tipping-point / partial-\(R^2\) sensitivity (sensitivity_report) (Cinelli and Hazlett 2020).
  5. Treat discovery graphs as sensitivity only (merge_discovery_sensitivity!), never as silent DAG replacement.

Package prose and DOIs: CausalTargeted methods Β· references. Worked Ξ΄-grid: Chapter 23 (Policy Evaluation).

37.3 Package layout

CausalDynamics.jl (identification and simulation):

CausalDynamics.jl/
β”œβ”€β”€ src/
β”‚   β”œβ”€β”€ CausalDynamics.jl
β”‚   β”œβ”€β”€ graphs/              # DAG ops, paths, CausalGraph, time_indexed
β”‚   β”œβ”€β”€ identification/      # Backdoor, frontdoor, IV, identify
β”‚   β”œβ”€β”€ scm/                 # GraphSCM, do(Β·), shared-U counterfactuals
β”‚   β”œβ”€β”€ cdm/                 # DiscreteTimeCDM trajectories
β”‚   β”œβ”€β”€ interventions/       # Abstract intervention types
β”‚   β”œβ”€β”€ integration/         # TMLE / RxInfer / discovery / SciML bridges
β”‚   └── utils/
β”œβ”€β”€ ext/                     # Weakdep extensions
β”œβ”€β”€ test/
β”œβ”€β”€ examples/
β”‚   β”œβ”€β”€ discrete_cdm.jl
β”‚   └── sciml_cdm_recipe.jl
β”œβ”€β”€ docs/
└── Project.toml

Continuous ODEs compose via optional SciML helpers (using OrdinaryDiffEq β†’ solve_cdm), not a solver hard-dependency in core. Deferred continuous-time / filter modules are listed on the package Scope page.

CausalTargeted.jl is the estimation companion (not nested inside CausalDynamics):

CausalTargeted.jl/
β”œβ”€β”€ src/
β”‚   β”œβ”€β”€ small_n.jl           # recommend_run_options, lean SL defaults
β”‚   β”œβ”€β”€ mtp_learners.jl      # Super Learner stacks
β”‚   β”œβ”€β”€ lmtp_grid.jl         # continuous MTP Ξ΄-grids
β”‚   β”œβ”€β”€ mediation_grid.jl    # interventional mediation under MTP
β”‚   β”œβ”€β”€ sequential_lmtp.jl
β”‚   β”œβ”€β”€ positivity.jl
β”‚   β”œβ”€β”€ sensitivity.jl
β”‚   └── …
β”œβ”€β”€ ext/                     # MLJ / EvoTrees / Flux weakdeps
β”œβ”€β”€ docs/src/
β”‚   β”œβ”€β”€ methods.md
β”‚   β”œβ”€β”€ references.md
β”‚   └── small_n.md
└── test/

Dual-stack synthetic recovery vs R lives in the book repo under scripts/synthetic_benchmark/ (treat single-draw wins as provisional).

Layout principles: separate graphs, identification, static SCM, and discrete CDMs; keep estimation out of CausalDynamics; prefer weakdep extensions over hard dependencies for optional backends.

37.4 End-to-end: identify β†’ LMTP

A thin application should identify on a DAG, then estimate with CausalTargeted, without baking nuisance fitting into CausalDynamics:

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 CausalTargeted Graphs DataFrames StableRNGs

g = DiGraph(3)
add_edge!(g, 1, 2)
add_edge!(g, 1, 3)
add_edge!(g, 2, 3)
id = identify(g, TotalEffectQuery(2, 3); node_names = [:W, :A, :Y])
println("Certificate: ", id.strategy, " adjustment = ", id.adjustment)

df, _ = simulate_linear_mtp(80; rng = StableRNG(30))
opts = recommend_run_options(nrow(df); engine = :lmtp)
grid = run_lmtp_grid(
    df, :A, :Y;
    baseline = [:W],
    deltas = [0.0, 0.5],
    folds = opts.folds,
    learners_outcome = DEFAULT_SL_LEARNERS,
    parallel = false,
    rng = StableRNG(31),
)
println("LMTP TE at Ξ΄ = ", grid.delta[end], ": ", round(grid.est[end]; digits = 3))
Certificate: backdoor adjustment = [:W]
LMTP TE at Ξ΄ = 0.5: 0.259

Static SCMs use the same stack’s GraphSCM / simulate_scm / do_intervention path (Chapter 4); discrete-time shared-\(U\) trajectories use DiscreteTimeCDM (Chapter 28).

37.5 Types and tests that matter here

The supported discrete-time type is approximately:

struct DiscreteTimeCDM{I, S, N} <: AbstractCDM
    variables::Vector{Symbol}
    initialise::I      # (rng) -> NamedTuple
    sample_noise::S    # (rng, state, t) -> NamedTuple
    step::N            # (state, t, noise, intervention) -> NamedTuple
end

Chapter 28 builds models with this type and calls simulate / counterfactual. Prefer package tests as the template: unit tests for a single mechanism, integration tests that intervene and replay shared noise, and synthetic recovery checks with known truth in CausalTargeted / the book benchmark scripts. Do not invent parallel APIs in application code that duplicate DoIntervention, DoSequence, or LMTP grids already exported by the packages.

37.6 Reproducibility

Pin Julia and major dependencies in Project.toml / Manifest.toml (this book targets Julia 1.12+ for CausalDynamics and CausalTargeted). Document which packages must be developed from path or Git URL versus General (DAGMakie is on General; see the package READMEs for the others). Examples in this book load dependencies via scripts/ensure_packages.jl.

37.7 Summary

The working stack separates identify (CausalDynamics), estimate (CausalTargeted), and display (DAGMakie). Application repositories should stay thin: data and concordance, not reimplemented estimators. Layout, weakdep extensions, and shared-\(U\) tests are the practices that keep that separation honest.

37.8 Further Reading

37.9 AgeSCM engine matrix

The AgeSCM reference application (MIRS mosquito age) illustrates combining packages without forcing one PPL to ingest thousands of wavenumbers.

Task Recommended engine Notes
Causal discovery (PC, OCE) Associations.jl Optional weakdep; bridge via prepare_from_discovery (Ch. 05b)
Causal graph & identification CausalDynamics.jl, DAGMakie Adjustment sets; RxInfer export
Discrete-time CDM trajectories CausalDynamics.jl (DiscreteTimeCDM) Ch. 28; shared-U CF
Continuous MTP / mediation CausalTargeted.jl LMTP Ξ΄-grids, positivity, sensitivity (Ch. 20–23)
Model specification GraphPPL.jl @model on αΊ‘, not raw Y
Scalable inference RxInfer.jl Variational message passing (AgeSCM: installed; wiring demo on αΊ‘)
Spectral encoding Flux.jl 1D CNN β†’ αΊ‘ ∈ ℝ^d
Domain shift IPTW (+ TMLE.jl Ch 20) Propensity of C; T2 mean OLS MAE 22.87 (full LOSO)
Small-n checks Turing.jl Optional NUTS
TipWhen to use RxInfer (VI) versus Turing (MCMC)

Both engines answer Bayesian questions; they differ in model class and cost.

Prefer RxInfer / GraphPPL (variational message passing) when:

  • the generative head is (or can be reduced to) a conjugate / factorised model after residualisation (e.g. CausalDynamics infer_backdoor_effect on demeaned \((y,x)\))
  • \(N\) is large and you need a fast posterior summary for a scalar effect \(\tau\) or a low-dimensional latent
  • you already sit on an AgeSCM-style pipeline: encoder β†’ αΊ‘ β†’ GraphPPL head

Prefer Turing (NUTS or other MCMC) when:

  • the likelihood is non-conjugate, hierarchical, or otherwise awkward for the current GraphPPL factors (abduction soft likelihoods, partial pooling, experimental-design comparisons)
  • you need a joint posterior over several parameters, or a small-\(n\) check against a VI mean
  • you want posterior predictive replicates drawn from full MCMC (Ch. 21)

Also in the stack (not a PPL choice): for latent trajectories, prefer filters / smoothers / StateSpaceDynamics Laplace before either PPL; use Turing only on short series or for parameter blocks once the path is handled.

Rule of thumb: conjugate tabular ATE and belief updates β†’ RxInfer; flexible or pedagogical continuous models β†’ Turing; long dynamical latent paths β†’ SSM tooling first. Local timings live in scripts/bayesian_benchmarks.jl.

Design rule: αΊ‘ is the PPL observation interface; TMLE-style propensity supports transport-weighted training. Resolve Graphs ≀1.13 with RxInfer 5.x (DataStructures 0.18) in AgeSCM Project.toml, see AgeSCM docs/RXINFER_DEPS.md.

37.9.1 Live wiring: CausalDynamics β†’ RxInfer

The book environment now ships GraphPPL and RxInfer. A confounded Gaussian table recovers \(\tau\approx 2\) in a fraction of a second:

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

Random.seed!(42)
n = 500
z = randn(n)
x = z .+ 0.4 .* randn(n)
y = 2.0 .* x .+ z .+ 0.3 .* randn(n)
data = DataFrame(Z = z, X = x, Y = y)
g = DiGraph(3)
add_edge!(g, 1, 2); add_edge!(g, 1, 3); add_edge!(g, 2, 3)
names = Dict(1 => :Z, 2 => :X, 3 => :Y)

t_w = @elapsed begin
    global wire = infer_backdoor_effect(g, data, 2, 3; node_names = names, iterations = 30)
end
println("Identifiable: ", wire.identifiable,
    "  confounders = ", wire.confounders,
    "  Ο„_mean β‰ˆ ", round(wire.Ο„_mean; digits = 3),
    "  (", round(t_w; digits = 3), " s)")
Identifiable: true  confounders = [:Z]  Ο„_mean β‰ˆ 2.001  (8.821 s)

AgeSCM replaces raw spectra by an encoder αΊ‘ before this head; the identification and VI steps stay the same. Timings: scripts/bayesian_benchmarks.jl.