7  Causal Discovery

Causal Discovery

Status: Draft

v0.1

7.1 From Identification to Discovery

ImportantCentripetal interlude (Observable → Structural)

Part I is otherwise depth-first: we assume structure, then ask what data can identify. This chapter reverses the arrow: infer candidate graphs from the Observable stratum. Treat discovered structure as a hypothesis / sensitivity input to identification—not as a silent replacement for a scientifically argued DAG (Chapter 29; CausalTargeted merge_discovery_sensitivity!).

Previous chapters assumed that the causal graph \(G\) was known. Identification theory (see Identification: When Can We Learn from Data?) determines when a causal question can be answered under that assumption. In practice, the graph may itself be uncertain and must be inferred from observational, and sometimes interventional, data. This chapter examines causal discovery (Spirtes et al. 2000; Pearl 2009).

The fundamental challenge is that data alone cannot distinguish all causal structures. Consider three variables \(X\), \(Y\), and \(Z\). The graphs \(X \rightarrow Y \rightarrow Z\), \(X \leftarrow Y \leftarrow Z\), and \(X \leftarrow Y \rightarrow Z\) may all imply the same set of conditional independencies under the Markov assumption. Such graphs form a Markov equivalence class: they are indistinguishable from observational data. Discovery algorithms typically recover this equivalence class rather than a unique DAG.

From the perspective of the three strata (see The Causal Hierarchy and Three Strata), discovery proposes candidates for the Structural stratum from patterns in the Observable stratum. A discovered graph is therefore a structural hypothesis, not a direct recovery of ontology. Only after scientific review and suitable assumptions does it support interventional or counterfactual reasoning in the CDM framework.

7.2 Constraint-Based Methods

Constraint-based methods infer causal structure by testing conditional independencies in the data. The key insight: if \(X \perp\!\!\!\perp Y \mid Z\) holds in the population, then no direct edge exists between \(X\) and \(Y\) in the true DAG (under faithfulness). By systematically testing such independencies, we can prune edges and orient others.

7.2.1 The PC Algorithm

The PC algorithm (Spirtes et al. 2000; JuliaDynamics 2024) is the canonical constraint-based method. It proceeds in two phases:

  1. Skeleton identification: Start with a complete undirected graph. For each pair \((X, Y)\), test \(X \perp\!\!\!\perp Y \mid S\) for conditioning sets \(S\) of increasing size. If independence holds for some \(S\), remove the edge \(X\), \(Y\).
  2. Orientation: Apply orientation rules (e.g., v-structures) to direct edges using the pattern of conditional independencies.

The algorithm relies on the faithfulness assumption: the conditional independencies in the data exactly correspond to d-separation in the true DAG. No “accidental” independencies exist beyond those implied by the graph structure.

Pseudocode (conceptual):

PC(D, α):
  1. Form complete undirected graph C on vertices V
  2. n ← 0
  3. repeat:
  4.   for each pair (X,Y) with edge in C:
  5.     for each S ⊆ Adj(X)\{Y} with |S| = n:
  6.       if X ⊥ Y | S (test at level α):
  7.         Remove edge X—Y, record S as SepSet(X,Y)
  8.         break
  9.   n ← n + 1
  10. until no changes
  11. Orient v-structures: X—Z—Y, X—/—Y ⇒ X→Z←Y
  12. Apply further orientation rules (Meek rules)
  13. return (partially directed) graph

Complexity: The PC algorithm has complexity \(O(p^d)\) where \(p\) is the number of variables and \(d\) is the maximum degree of the true graph. Sparse graphs admit efficient discovery.

The following example uses Associations.jl for PC on a simulated observational cohort with nutrition confounding treatment and worm burden (nutrition → treatment → worm_burden, nutrition → worm_burden), then passes the candidate graph to CausalDynamics for backdoor identification (JuliaDynamics 2024; Pedersen and Babayan 2011).

WarningCPDAG vs DAG

PC returns a partially directed graph (CPDAG): some edges may remain bidirectional within a Markov equivalence class. prepare_from_discovery accepts the candidate structure; use complete=true to break bidirectional edges before backdoor adjustment when a fully oriented DAG is required.

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

include(joinpath(dirname(Base.active_project()), "scripts", "biological_showcases.jl"))

cohort_data = confounded_cohort_dgp(; n = DEFAULT_COHORT_N, rng = StableRNG(DEFAULT_COHORT_SEED))
2500×3 DataFrame
2475 rows omitted
Row nutrition treatment worm_burden
Float64 Float64 Float64
1 -0.670252 -0.575533 -0.769374
2 0.447122 0.633522 0.483733
3 1.37363 1.24678 2.30408
4 1.30954 0.757185 1.06823
5 0.12607 0.878857 0.668735
6 0.683948 0.438237 0.72563
7 -1.0192 -0.767335 -1.0106
8 -0.793513 -0.537246 -1.11888
9 1.77472 1.5511 1.81449
10 1.29735 1.09219 1.52236
11 -1.64385 -1.37354 -1.51247
12 0.794439 0.267494 0.410471
13 -1.30967 -1.25093 -1.51639
⋮ ⋮ ⋮ ⋮
2489 3.57935 2.59069 3.40992
2490 1.48404 0.918325 1.09554
2491 1.58608 1.58189 2.35057
2492 -0.174721 -0.337562 -0.0998417
2493 0.799415 0.90453 0.892252
2494 0.22084 0.0534283 0.0372172
2495 -0.00238338 -0.12911 -0.15881
2496 0.276069 -0.723027 -0.522364
2497 0.146124 0.308617 0.552386
2498 3.05265 2.77732 3.80994
2499 0.426271 -0.0166701 0.253988
2500 -1.02834 -0.580304 -0.584273 
include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using CausalDynamics Associations

ĝ = infer_pc_graph(cohort_data, [:nutrition, :treatment, :worm_burden]; verbose = false)
confounders, identifiable = prepare_from_discovery(ĝ, :treatment, :worm_burden; complete = true)
@assert identifiable
@assert :nutrition in confounders
pc_names = [:nutrition, :treatment, :worm_burden]
(confounders, identifiable)
([:nutrition], true)

PC on the confounded cohort. Left: undirected skeleton of the discovered graph (dagplot_skeleton). Right: identifiable fork for treatment → worm_burden with nutrition highlighted as the recovered confounder.
NoteDiscovery is not identification

Recovering a graph from data yields a hypothesis under faithfulness and sufficiency. Identification (Ch. 5) then asks whether a causal effect is computable from that structure. Conceptually, discovery and identification are separate steps; in these examples, Associations.jl proposes structure and CausalDynamics issues identification certificates. In production pipelines, treat a discovered adjustment set as sensitivity relative to a user DAG, helpers such as CausalTargeted merge_discovery_sensitivity! record disagreement and must not silently replace the analyst’s graph (Reporting Standards).

In observational cohorts (anthelmintic, vaccine, or otherwise), host or contextual condition often confounds who receives treatment and subsequent outcomes, the PC step above mirrors that pattern with generic column names. For longitudinal markers and load, lagged coupling is a natural target for OCE and temporal backdoor adjustment.

Conceptual PC skeleton (pseudocode) 
PC(D, α):
  1. Form complete undirected graph C on vertices V
  2. n ← 0
  3. repeat:
  4.   for each pair (X,Y) with edge in C:
  5.     for each S ⊆ Adj(X)\{Y} with |S| = n:
  6.       if X ⊥ Y | S (test at level α):
  7.         Remove edge X—Y, record S as SepSet(X,Y)
  8.         break
  9.   n ← n + 1
  10. until no changes
  11. Orient v-structures: X—Z—Y, X—/—Y ⇒ X→Z←Y
  12. Apply further orientation rules (Meek rules)
  13. return (partially directed) graph

7.2.2 The FCI Algorithm

When latent confounders may exist, the PC algorithm can produce incorrect conclusions. The FCI (Fast Causal Inference) algorithm (Spirtes et al. 2000; Zhang 2008) extends PC to handle latent variables and selection bias.

FCI outputs a Partial Ancestral Graph (PAG) rather than a DAG. A PAG represents the Markov equivalence class of the underlying DAG with latent variables. Edge marks include:

  • \(-\) (solid): definite adjacency
  • \(o\) (circle): possibly arrowhead or tail (uncertain orientation)
  • \(>\) (arrowhead): definite arrowhead

When to use FCI vs PC:

  • PC: Use when you assume causal sufficiency, no unobserved common causes. Appropriate for carefully designed experiments or when confounders are measured.
  • FCI: Use when latent confounders are plausible. FCI correctly represents uncertainty about orientation and can indicate possible confounding.

Interpreting a PAG requires care: a circle (\(o\)) indicates that the algorithm could not determine whether the edge is directed or bidirected from the data. Such ambiguity reflects genuine epistemic limits, observational data cannot resolve certain structural questions without further assumptions or interventions.

7.3 Score-Based Methods

Score-based methods treat structure learning as optimisation: search over DAGs to maximise a score that balances fit to data against model complexity.

The general form of the score is:

\[ \text{score}(G, D) = \log \mathcal{L}(D \mid G, \hat{\theta}) - \text{penalty}(|G|) \]

where \(\hat{\theta}\) are the maximum-likelihood parameters for \(G\) given data \(D\), and the penalty discourages overly complex graphs.

7.3.1 BIC and AIC

Common scores include:

  • BIC (Bayesian Information Criterion): \(\text{BIC} = -2 \log \mathcal{L} + k \log n\), where \(k\) is the number of parameters and \(n\) is sample size. BIC is consistent for structure selection under regularity conditions.
  • AIC: \(\text{AIC} = -2 \log \mathcal{L} + 2k\). AIC tends to select larger models; BIC penalises complexity more strongly.

7.3.2 Greedy Equivalence Search (GES)

GES (Chickering 2002) operates on the space of equivalence classes (represented as completed partially directed acyclic graphs, CPDAGs) rather than individual DAGs. It performs a two-phase greedy search:

  1. Forward phase: Start with empty graph. Repeatedly add the edge that most improves the score.
  2. Backward phase: Repeatedly remove the edge that most improves the score.

GES is score-equivalent: it assigns the same score to Markov-equivalent DAGs. Under faithfulness and with a consistent score (e.g., BIC), GES is consistent for the true equivalence class.

7.3.3 Bayesian Approach

A fully Bayesian approach places a posterior over DAGs:

\[ P(G \mid D) \propto P(D \mid G) \, P(G) \]

where \(P(D \mid G)\) is the marginal likelihood (integrated over parameters) and \(P(G)\) is a prior over structures. Common priors favour sparser graphs (e.g., uniform over edge count, or a Beta-Binomial prior). MCMC or exact methods can approximate the posterior. This yields not just a point estimate but uncertainty over structures, useful when multiple graphs fit the data similarly.

Comparison with constraint-based methods:

Aspect Constraint-based (PC) Score-based (GES)
Assumptions Faithfulness, causal sufficiency Faithfulness, parametric model
Output Equivalence class (CPDAG) Equivalence class (CPDAG)
Robustness Sensitive to CI test errors Sensitive to model misspecification
Scalability \(O(p^d)\) Depends on search; often faster in practice

7.4 Time Series Causal Discovery

Causal discovery in time series must account for temporal ordering and autocorrelation. The past can cause the present, but not vice versa, this constraint reduces the space of possible structures.

7.4.1 Granger Causality

Granger causality (Granger 1969) is a predictive notion: \(X\) Granger-causes \(Y\) if past values of \(X\) improve prediction of \(Y\) beyond past values of \(Y\) alone.

Mathematical formulation (bivariate case): Let \(Y_t\) be predicted by:

\[ Y_t = \sum_{k=1}^{K} \alpha_k Y_{t-k} + \sum_{k=1}^{K} \beta_k X_{t-k} \]

\(X\) Granger-causes \(Y\) if any \(\beta_k \neq 0\) (reject the null that all \(\beta_k = 0\)).

Limitations:

  • Confounders: A common cause \(Z\) of \(X\) and \(Y\) can induce apparent Granger causality.
  • Nonlinearity: Granger causality is linear; nonlinear dependencies may be missed.
  • Indirect causation: Granger detects temporal precedence, not necessarily direct causation.

7.4.2 PCMCI

PCMCI (Peter and Clark Momentary Conditional Independence) (Runge et al. 2019) combines the PC algorithm with time-lagged dependencies. It:

  1. Builds a superset of parents using a variant of PC with lagged variables.
  2. Applies MCI (Momentary Conditional Independence) tests to remove spurious links, conditioning on the full past.

PCMCI properly handles autocorrelation: conditioning on past values of all variables avoids spurious associations from shared history. The key assumption is causal sufficiency in the time-lagged variable set, no unobserved common causes of the lagged variables.

7.4.3 Optimal causation entropy (OCE)

Optimal causation entropy (OCE) (Sun et al. 2015) selects lagged parents for each series via pairwise and conditional association tests. In Julia, Associations.jl implements OCE and returns parent sets that convert to CausalDynamics TemporalDAGSpec for unrolled identification (see Ch. 28) (JuliaDynamics 2024).

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

include(joinpath(dirname(Base.active_project()), "scripts", "biological_showcases.jl"))

ts = immune_parasite_ts_dgp(; T = DEFAULT_TS_LENGTH, rng = StableRNG(DEFAULT_TS_SEED))
spec = infer_oce_temporal_spec(ts.series, ts.variables; verbose = false)
u = unroll_temporal_dag(spec, 5)
adj_nodes = temporal_backdoor_adjustment_nodes(u, :immune, 1, :parasite_load, 2)
(spec.variables, length(spec.edges), adj_nodes)
([:immune, :parasite_load], 4, Set([(:parasite_load, 0), (:immune, 0)]))

OCE lag parents unrolled over five occasions (dagplot_temporal). Columns are time; rows are immune and parasite load. Feedback and autoregressive lags appear as edges across neighbouring columns.
NoteOCE lag convention

Associations stores embedding lags in parents_τs (negative indices). CausalDynamics maps these to LaggedEdge lags with lag = abs(τ), so τ = -1 becomes \(X_{t-1} \rightarrow Y_t\) under unroll_temporal_dag.

7.4.4 Convergent Cross Mapping (CCM)

CCM (Sugihara et al. 2012) is designed for nonlinear dynamical systems. It exploits Takens’ embedding theorem: the attractor of a dynamical system can be reconstructed from a single observed time series. CCM tests whether \(Y\)’s attractor can be reconstructed from \(X\)’s time series, if so, \(X\) and \(Y\) are dynamically coupled.

7.4.5 Interventional Dynamical Causality (IntDC)

Granger causality, transfer entropy, and CCM quantify what Shi and colleagues call constructive dynamical causality (ConDC) (Shi et al. 2026): dependence of a target’s evolution on a putative cause in the observed dynamics (prediction or reconstructability). That is distinct from interventional dynamical causality (IntDC), which asks how the target would change under an intervention on the cause.

When ethical or practical constraints rule out real perturbations, IntDC must be inferred from observational series alone. Shi et al. propose Interventional Embedding Entropy (IEE) (Shi et al. 2026): an information-flow measure in delay-embedding space that approximates interventional response from local neighbourhood geometry, without requiring a fitted dynamical model or experimental interventions. Relative to ConDC indices, IEE is intended to rank edge importance and to reduce confounding-driven spurious links in multivariate reconstructions. Applications in the source paper include connectome estimation, epidemic transmission networks, and circadian regulatory structure. In the book’s pipeline, IEE sits with CCM and transfer entropy as an embedding-space screen; candidate edges still need identification and estimation under an explicit SCM or CDM (Chapters 5 and 28). CausalDynamics ships a Julia port of the reference algorithm (smsxiaomayi/IEE); with Associations loaded, mi = :auto uses KSG1 mutual information, while mi = :reference keeps the MATLAB-faithful estimator used in concordance tests.

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

# Coupled logistic (x → y): IEE ranks the true direction, then feeds a TemporalDAGSpec
rng = Random.Xoshiro(1)
T = 400
x = zeros(T); y = zeros(T)
x[1] = 0.1; y[1] = 0.2
β = 0.15
for t in 1:(T - 1)
    x[t + 1] = 3.7 * x[t] * (1 - x[t]) + 0.01 * randn(rng)
    y[t + 1] = 3.7 * y[t] * (1 - (1 - β) * y[t] - β * x[t]) + 0.01 * randn(rng)
end
x = x[101:end]; y = y[101:end]

iee_xy = interventional_embedding_entropy(x, y; p = 2, k = 2, n_delta = 8, mi = :reference)
iee_yx = interventional_embedding_entropy(y, x; p = 2, k = 2, n_delta = 8, mi = :reference)
scores = iee_score_matrix([x, y]; p = 2, k = 2, n_delta = 8, mi = :reference)
thr = 0.5 * (scores[1, 2] + scores[2, 1])
spec = iee_to_temporal_spec(scores, [:x, :y]; threshold = thr, lag = 1)
u = unroll_temporal_dag(spec, 4)
(iee_xy, iee_yx, length(spec.edges), nv(u.graph))
(0.18462368168744026, 0.027359979226589986, 1, 10)

IEE on coupled logistics (true x \rightarrow y). Left: series. Centre: pairwise IEE scores (heatmap; diagonal blank). Right: thresholded temporal DAG via dagplot_temporal (occasions \(t = 0:T\) with \(T = 4\)).

When to use:

  • Granger: Linear systems, quick screening, well-established interpretation. Suitable when relationships are approximately linear and confounders are controlled.
  • PCMCI: Multivariate time series, proper handling of autocorrelation, constraint-based approach. Preferred when multiple variables are observed and lagged dependencies matter (see tigramite for a reference Python implementation).
  • OCE: Julia-native lag parent selection via Associations.jl; composes with time-indexed identification in the book’s examples.
  • CCM: Nonlinear dynamics, potentially low-dimensional attractors, when linear methods may fail. Particularly useful in ecology and climate science where systems exhibit nonlinear coupling.
  • IEE / IntDC: When the scientific target is interventional ranking from observational series (non-intervention systems), and ConDC screens risk confounding or fail to order effect strength (Shi et al. 2026). CausalDynamics exports interventional_embedding_entropy / infer_iee_temporal_spec as a bridge into time-indexed identification (mi = :auto uses Associations KSG1 when loaded; mi = :reference matches the MATLAB port).
  • ODE parents across environments: When multiple experimental environments are available and the target is an ODE mechanism for one response that should remain stable across those environments (Pfister et al. 2019; Peters et al. 2022). Use infer_ode_parents / ode_parent_ranking_to_continuous_spec (DataInterpolations for spline derivatives). Full constrained-spline scoring from the CausalKinetiX R reference remains future work.

Note that Granger causality and structural causation are distinct concepts: \(X\) may Granger-cause \(Y\) without being a direct cause (e.g., through a mediating variable), and \(X\) may cause \(Y\) without Granger-causing it (e.g., in nonlinear systems where the relationship is not captured by linear prediction). ConDC versus IntDC sharpens the same point for dynamical discovery: predictive coupling is not automatically an interventional effect.

7.5 Discovery Under Interventions

Interventional data provides strictly more information than observational data (Eberhardt et al. 2007). An intervention on \(X\) breaks incoming edges to \(X\), altering the observational distribution. By comparing distributions under different interventions, we can distinguish graphs within the same Markov equivalence class. When interventions are unavailable, IntDC-style criteria such as IEE attempt to recover interventional rankings from observational embeddings (Shi et al. 2026), complementary to (not a substitute for) experimental design.

Experimental design for structure learning:

  • Adaptive interventions: Choose which variable to intervene on next based on current uncertainty about the graph.
  • Minimum intervention sets: Find the smallest set of interventions needed to identify the true DAG.

Active learning of causal structure iteratively selects interventions to maximise information gain, reducing the equivalence class until the true structure is identified. In the best case, \(O(\log p)\) interventions suffice to identify a DAG on \(p\) variables; in the worst case, \(p-1\) interventions may be needed. The structure of the equivalence class determines the minimum number required.

7.6 Assumptions and Limitations

Causal discovery rests on several assumptions that may not hold in practice:

7.6.1 Core Assumptions

  • Faithfulness: All conditional independencies in the population arise from d-separation in the true DAG. Violations (e.g., near-cancellation of paths, or fine-tuned parameters that create “accidental” independencies) can mislead algorithms into removing true edges or failing to orient correctly.
  • Causal sufficiency: No unobserved common causes. Violations require FCI or similar methods that output PAGs.
  • Acyclicity: The true causal structure is a DAG. Feedback loops require extensions (e.g., summary graphs for stationary processes, or dynamic Bayesian networks with time-indexed variables).

7.6.2 What Can and Cannot Be Learned

From observational data alone:

  • We can typically recover the Markov equivalence class (CPDAG or PAG).
  • We cannot uniquely determine edge orientations that are equivalent under the Markov property.
  • We cannot distinguish causation from correlation without further assumptions (e.g., faithfulness, parametric restrictions).

7.6.3 Sample Complexity

Discovery algorithms require sufficient data for reliable conditional independence tests or score estimation. Sample complexity depends on:

  • Graph sparsity
  • Effect sizes (weaker effects require more data)
  • The number of variables \(p\)

7.6.4 Connection to Identifiability

Identification (Identification: When Can We Learn from Data?) asks: given a graph, can we identify a causal effect? Discovery asks: can we identify the graph itself? The two are complementary: discovery provides the structural input for identification; identification tells us what we can learn once we have (or assume) a structure. In practice, analysts often iterate: use discovery to propose a graph, check identification for the effects of interest, then refine the graph or collect interventional data if identification fails.

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

# Observable → Associations (structure) → CausalDynamics (identification) → estimation
confounders_pc, ok_pc = discover_and_prepare(
    cohort_data, :treatment, :worm_burden;
    method = :pc,
    names = [:nutrition, :treatment, :worm_burden],
    complete = true,
    verbose = false,
)
# Downstream: TMLE / RxInfer / g_computation on the adjustment set (Ch. 5, Ch. 28)
(confounders_pc, ok_pc)
([:nutrition], true)

7.7 Connection to the CDM Framework

Causal discovery is the bridge from the Observable stratum to the Structural stratum. Discovered structures feed the CDM (Causal Dynamical Model) framework:

  1. Discovery yields a candidate DAG (or equivalence class) from data.
  2. Modelling instantiates the structure with parametric or nonparametric dynamical laws.
  3. Validation checks consistency with new data and domain knowledge.
  4. Refinement iterates: discovery → model → validate → discover.

The discovered structure enables interventional reasoning (what if we set \(X = x\)?) and counterfactual reasoning (what alternative trajectory follows under a declared cross-world coupling?). Without a structural hypothesis, these questions cannot be addressed. Discovery provides that hypothesis, with appropriate epistemic humility about equivalence classes and assumptions; the process vocabulary is a secondary interpretation (Ch. 9).

In the process-philosophy framing of this book, discovery proposes the structure through which occasions may be related, using observable traces as evidence. The resulting graph is more than a predictive summary when its causal interpretation is justified, but the discovery output alone does not establish the relation type, its mechanism, or the ontology of the entities represented.

Discovery is thus an inverse problem: we observe the outer stratum (associations, conditional independencies) and infer the inner Structural stratum. Like all inverse problems, it is ill-posed without assumptions. Faithfulness, causal sufficiency, and acyclicity provide the regularisation that makes the problem tractable. When these assumptions hold approximately, discovery algorithms offer a principled route from data to structure, and from structure to the interventional and counterfactual reasoning that the CDM framework enables.