35 Software Engineering Patterns for Causal-Dynamical Modelling

Status: Draft

v0.4

35.1 Introduction

Part IV needs a software practice companion. This chapter describes how the book’s Julia stack is organised, how a thin application connects identification to estimation, and how tests and dependency pins keep work reproducible. It is not a substitute for the conceptual chapters or the package manuals. The principles transfer to other languages; the concrete layout uses Julia. Together with reporting standards, it completes identify → estimate → display after the CDM synthesis.

35.2 Working stack

In these pages the working stack is four packages:

CausalDynamics.jl — graphs, identify certificates, SCMs / CDMs
CausalTargeted.jl — cross-fitted LMTP / MTP total-effect grids, sequential and thin survival LMTP, factor recodes, Super Learner profiles, positivity / sensitivity
CausalMediation.jl — TE / NDE / NIE and related mediated contrasts (run_mediation, moc)
DAGMakie.jl — the same graph as a figure

Exported names for identification and estimation match the literature (identify, do_intervention, run_lmtp_grid, run_mediation, dagplot, …). Temporal identity may use the closed lexicon in Table 8 (occasion, enduring, constitution, constitutive, influence, onset, replacement, deployment). The philosophical glossary belongs in this book (Chapter 9; glossary).

Division of labour:

Layer	Owns	Does not own
Identification / CDM simulation (CausalDynamics)	Graphs, `identify`, temporal unrolling, hierarchical nesting (`RandomEffectSpec`, `unroll_hierarchical_dag`), `GraphSCM` / `DiscreteTimeCDM`, `simulate_panel`	Cross-fitted nuisance grids
Estimation — LMTP (CausalTargeted)	Non-mediated LMTP / MTP grids, sequential and survival LMTP, factor recodes, positivity, sensitivity, `choose_policy`, transport weights	Mediated EIF / NDE–NIE
Estimation — mediation (CausalMediation)	Mediation grids, `moc`, RI / natural / RT when needed	Cohort loaders; non-mediated LMTP
Display (DAGMakie)	Adjustment / SWIG / temporal layouts	Estimators
Application repos	Data, registries, concordance, domain labels for transport	Re-implementing EIF estimators

On General, pin CausalDynamics ≥ 0.4 and CausalTargeted ≥ 0.3.10 (General#165016; 0.3.7–0.3.9 skipped). CausalMediation 0.1.1 (factor-A interventional mediation) registers next from the package tip once AutoMerge accepts it. The book develops the packages/ checkouts. Soft façades on CausalTargeted still forward to CausalMediation when that package is loaded.

When extending the book or the packages, keep Pearl/SciML identification and estimation names stable, allow the closed temporal lexicon on existing types, keep Whitehead glossary in book prose, and update Table 8 if you introduce new central mappings. For manuals and changelogs, use the repositories, not this chapter as API documentation.

Small-n checklist (CausalTargeted)

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

Start from recommend_run_options(n; engine, n_mediators) (lean Super Learner when \(n < 80\); parallel=false by default).
Pass mixed-type baselines as-is (String / Bool / categorical are encoded with a fold-stable CovariateSchema in LMTP, sequential, survival, and mediation grids).
Inspect positivity / support (positivity_report or grid positivity=true) (Petersen et al. 2012).
For mediation, sweep nested MC (mediation_n_mc_sweep) until signs and SEs stabilise (Liu et al. 2024).
Report tipping-point / partial-\(R^2\) sensitivity (sensitivity_report) (Cinelli and Hazlett 2020).
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).

35.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, simulate_panel
│   ├── 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 LMTP estimation companion (not nested inside CausalDynamics): cross-fitted total-effect grids, sequential / survival LMTP, factor recodes (DiscreteTreatmentPolicy / run_discrete_lmtp and sequential policies), positivity, and run_gcomp with a refitting bootstrap (or n_boot=0 for an influence-function SE). Incomplete outcomes and covariates use Observable policies: handle_missing (:drop, :ipcw, :impute, :ipcw_impute) on LMTP, g-comp, sequential, and survival paths, with run metadata for strategy and PCH rung. Structural MAR/MNAR claims attach via CausalDynamics MissingnessSpec / identify(...; missingness=); opt-in posterior draws use impute_posterior and run_lmtp_grid(...; imputation=). Survival censoring IPCW is not the same object as MAR missing terminal \(S_T\). Book overview: Observational Methods. Package methods (causal–statistical voice): CausalDynamics Missingness, CausalTargeted Missingness.

CausalTargeted.jl/
├── src/
│   ├── small_n.jl           # recommend_run_options, lean SL defaults
│   ├── mtp_learners.jl      # Super Learner stacks
│   ├── lmtp_grid.jl         # continuous MTP δ-grids (total effects)
│   ├── sequential_lmtp.jl
│   ├── discrete_lmtp.jl     # factor recodes (T=1)
│   ├── query_bridge.jl      # estimand_from_query
│   ├── survival_lmtp.jl     # discrete-time event-time / survival LMTP
│   ├── positivity.jl
│   ├── sensitivity.jl
│   └── …                    # soft mediation façades when CausalMediation loads
├── ext/                     # MLJ / EvoTrees / Flux / CausalMediation weakdeps
├── docs/src/
│   ├── methods.md
│   ├── references.md
│   └── small_n.md
└── test/

CausalMediation.jl owns mediated contrasts:

CausalMediation.jl/
├── src/
│   ├── identify_bridge.jl   # plan_mediation, spec_from_identification
│   ├── mediation_grid.jl    # run_mediation / run_mediation_grid
│   └── …
├── docs/
└── Project.toml

Dual-stack synthetic recovery vs R lives in the book repo under scripts/synthetic_benchmark/ (treat single-draw wins as provisional). Broader functionality, edge-case, and performance stress (mixed-type baselines, panels, small public tables) lives under scripts/stress_harness/ with the dataset catalog in data/catalog.toml. Methods write-up and dataset notebook: CausalTargeted stress validation (STRESS.md).

Layout principles: separate graphs, identification, static SCM, and discrete CDMs; keep estimation out of CausalDynamics; split LMTP (CausalTargeted) from mediation (CausalMediation); prefer weakdep extensions over hard dependencies for optional backends.

35.4 End-to-end: identify → LMTP

A thin application should identify on a DAG, then estimate with CausalTargeted (or CausalMediation for mediated contrasts; Ch. 20), without baking nuisance fitting into CausalDynamics. For discrete-time CDMs, prefer simulate_panel → sequential_spec_from_identification → execute_estimand (Ch. 28), with policies for factor \(A_t\); for event-free probability under MTP use SurvivalPolicy / run_survival_lmtp (same chapter; Ch. 23). The point-treatment path below remains the tabular default:

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)

rng = StableRNG(30)
df, truth_mix = simulate_mixed_baseline_mtp(80; rng = rng)
opts = recommend_run_options(nrow(df); engine = :lmtp)
grid = run_lmtp_grid(
    df, :A, :Y;
    baseline = truth_mix.baseline,
    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),
    " (mixed baseline encoded internally)")

Certificate: backdoor adjustment = [:W]
LMTP TE at δ = 0.5: 0.343 (mixed baseline encoded internally)

Static SCMs use the same stack’s GraphSCM / simulate_scm / do_intervention path (Chapter 4); discrete-time shared-\(U\) trajectories use DiscreteTimeCDM (Chapter 28). Book examples that need the host-immunity spine should load scripts/biological_showcases.jl (host_immunity_cdm, simulate_host_immunity_scenario, host_immunity_teaser_dag) rather than re-implementing the assignments inline.

35.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 / CausalMediation / the book benchmark scripts. Do not invent parallel APIs in application code that duplicate DoIntervention, DoSequence, LMTP grids, or mediation runners already exported by the packages.

35.6 Reproducibility

Pin Julia and major dependencies in Project.toml / Manifest.toml (this book targets Julia 1.13). On General: CausalDynamics 0.4+, CausalTargeted 0.3.7+ (fold-stable mixed-type baselines), CausalMediation 0.1.0+, DAGMakie (figures). Path Pkg.develop remains fine for book work before a release lands. Examples load dependencies via scripts/ensure_packages.jl.

35.7 Summary

The working stack separates identify (CausalDynamics), estimate LMTP (CausalTargeted) or mediation (CausalMediation), and display (DAGMakie for graphs; plot_mtp_curve for δ-indexed effect curves). 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.

35.8 Further Reading

CausalDynamics documentation (scope, API, examples)
CausalTargeted documentation
CausalMediation documentation
DAGMakie documentation
Concept Reference Table 8
Julia package development and performance guides (language-level detail)

35.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`, `simulate_panel`)	Ch. 28; shared-`U` CF; sequential LMTP hand-off
Nested units / clusters	CausalDynamics.jl (`RandomEffectSpec`, `unroll_hierarchical_dag`); CausalTargeted.jl (`cluster=` on MSM)	Generative nest vs cluster-robust IF; LMM in apps only
Latent / continuous / transport	CausalDynamics + CausalTargeted	`ObservationBridge`; continuous `g_computation`; `TransportQuery` / `choose_policy`
Continuous MTP (total effect)	CausalTargeted.jl	LMTP δ-grids, positivity, sensitivity (Ch. 23)
Event-time / survival LMTP	CausalTargeted.jl	`SurvivalPolicy` / `run_survival_lmtp` (Ch. 23, 28); competing risks deferred
Mediation (NDE/NIE, `moc`)	CausalMediation.jl	Certificate → `run_mediation` (Ch. 20)
Spectral / image encoding	CausalDynamics `RepresentationSpec` (+ Lux or Flux encoder callable)	Codes ẑ are DAG/estimation nodes; CNN layers are not. AgeSCM may still train a Flux 1D CNN and pass `encode` into the bridge
Model specification	GraphPPL.jl	`@model` on ẑ, not raw Y
Scalable inference	RxInfer.jl	Variational message passing (AgeSCM: installed; wiring demo on ẑ)
Graph-constrained deep `f_i`	CausalDynamics Lux `MechanismLibrary`	ODE residuals / static / additive generative L3; see SciML recipes
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
Deep SCM stress	CausalDynamics + CausalTargeted Quarto	`docs/stress/deep_scm_*.qmd` (not the book)

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

Computational troubleshooting: if VI and NUTS disagree on a scalar \(\tau\), or NUTS shows unstable geometry on a small check, treat that as a fitting problem to diagnose (approximation class, scaling, multimodality) before expanding scientific claims (Gelman et al. 2026). Prefer recovering a known contrast from simulated CDM/SCM data first (Ch. 5, 21).

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.

35.9.1 Representation bridge in the AgeSCM pattern

AgeSCM trains a spectral encoder and treats ẑ as the observation interface. In the owned stack the same contract is RepresentationSpec / encode_to_panel: the CNN (or linear map) is the encoder callable; the code columns are what identification, mediation, and generative L3 see. The chunk below uses a matched filter for speed; swap in a Lux or Flux encode without changing the DAG.

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

rng = Random.Xoshiro(42)
n, p = 80, 16
grid = range(0, 1; length = p)
bump = exp.(-((collect(grid) .- 0.45) ./ 0.1) .^ 2)
A = Float64.(rand(rng, n) .< 0.5)
Mstar = 1.1 .* A .+ 0.3 .* randn(rng, n)
S = randn(rng, n, p) .* 0.1
for i in 1:n
    S[i, :] .+= Mstar[i] .* bump
end
wide0 = DataFrame(A = A, Y = 1.6 .* Mstar .+ 0.2 .* randn(rng, n))
spec = RepresentationSpec(:spectrum, [:z1], Sm -> hcat(Sm * bump); role = :definitional)
wide = encode_to_panel(wide0, S, spec)
println(representation_certificate(spec))
println("cor(z1, M*) ≈ ", round(cor(wide.z1, Mstar); digits = 3),
    "; columns for estimation: ", names(wide))

(role = :definitional, source = :spectrum, code_names = [:z1], encode_type = "Main.Notebook.var\"#10#11\"", encode_type_hash = 0xd875467050136982, n_codes = 1)
cor(z1, M*) ≈ 0.992; columns for estimation: ["A", "Y", "z1"]

Mediation and LMTP then run on wide with :z1 as mediator or baseline feature; see CausalMediation methods and CausalTargeted Deep SCM estimation stress. Graph-constrained f_i and additive L3 on those codes are Chapter 28’s Deep SCM chunk.

35.9.2 CDM and estimation heads as software

Package layout (above) is half the story. Generative CDMs, PPL heads, and LMTP pipelines behave like software under test (Gelman et al. 2026):

Oracle recovery: simulate from a known SCM/CDM, run identification + estimation, assert \(\hat\tau\) (or learned dynamics) matches truth within tolerance.
Golden contrasts: keep a small fixed synthetic table (as in the RxInfer wiring chunk below) in CI or stress notebooks so regressions surface early.
Out-of-domain splits: treat LOSO / temporal hold-outs as regression tests when transport is claimed (Ch. 21).

Owned-package stress notebooks under CausalTargeted (and CausalDynamics Deep SCM stress) already follow this pattern; book chapters show the scientific side of the same checks.

35.9.3 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  (7.963 s)

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