Previous chapters assumed the causal graph \(G\) was known. Identification theory (see Identification: When Can We Learn from Data?) tells us when we can answer causal questions given a known structure. But in practice, we often do not know the graph, we must learn it from data. This chapter addresses causal discovery: the problem of inferring causal structure from observational (and sometimes interventional) data (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 aims to recover the Structural stratum (the DAG encoding prehensive relations) from patterns in the Observable stratum. The discovered structure then enables interventional and counterfactual reasoning in the CDM framework.
8.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.
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\).
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.
Precompiling packages...
5425.8 ms ✓ CausalDynamics → CausalDynamicsDAGMakieExt
6000.4 ms ✓ DAGMakie → DAGMakieCausalInferenceExt
2 dependencies successfully precompiled in 7 seconds. 333 already precompiled.
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
8.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.
8.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.
where \(\hat{\theta}\) are the maximum-likelihood parameters for \(G\) given data \(D\), and the penalty discourages overly complex graphs.
8.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.
8.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:
Forward phase: Start with empty graph. Repeatedly add the edge that most improves the score.
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.
8.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
8.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.
8.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:
\(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.
8.4.2 PCMCI
PCMCI (Peter and Clark Momentary Conditional Independence) (Runge et al. 2019) combines the PC algorithm with time-lagged dependencies. It:
Builds a superset of parents using a variant of PC with lagged variables.
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.
8.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).
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.
8.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.
8.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 TemporalDAGSpecrng =Random.Xoshiro(1)T =400x =zeros(T); y =zeros(T)x[1] =0.1; y[1] =0.2β =0.15for t in1:(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)endx = 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, 8)
IEE on coupled logistics (true x \rightarrow y). Left: series. Centre: pairwise IEE scores (heatmap; diagonal blank). Right: thresholded temporal DAG unrolled over four occasions.
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.
8.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.
8.6 Assumptions and Limitations
Causal discovery rests on several assumptions that may not hold in practice:
8.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).
8.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).
8.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\)
8.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.
Causal discovery is the bridge from the Observable stratum to the Structural stratum. Discovered structures feed the CDM (Causal Dynamical Model) framework:
Discovery yields a candidate DAG (or equivalence class) from data.
Modelling instantiates the structure with parametric or nonparametric dynamical laws.
Validation checks consistency with new data and domain knowledge.
Refinement iterates: discovery → model → validate → discover.
The discovered structure enables interventional reasoning (what if we set \(X = x\)?) and counterfactual reasoning (what alternative concrescence would follow for one organism with fixed creative advance\(\mathbf{u}\)? See Ch. 9). Without a structural hypothesis, these questions cannot be addressed. Discovery provides that hypothesis, with appropriate epistemic humility about equivalence classes and assumptions.
In the process philosophy framing of this book, discovery aims to recover the prehensive structure (the pattern of how occasions prehend one another) from its observable traces. The structure is not merely a convenient summary; it encodes the causal architecture that constrains what can happen under intervention and counterfactual supposition.
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.
Chickering, David Maxwell. 2002. “Optimal Structure Identification with Greedy Search.”Journal of Machine Learning Research 3: 507–54.
Eberhardt, Frederick, Clark Glymour, and Richard Scheines. 2007. “On the Number of Experiments Required to Identify the Causal Structure of a System.”Journal of Machine Learning Research 8: 2651–98.
Granger, Clive W. J. 1969. “Investigating Causal Relations by Econometric Models and Cross-Spectral Methods.”Econometrica 37 (3): 424–38.
Peters, Jonas, Stefan Bauer, and Niklas Pfister. 2022. “Causal Models for Dynamical Systems.” In Probabilistic and Causal Inference: The Works of Judea Pearl, edited by Hector Geffner, Rina Dechter, and Joseph Y. Halpern. Association for Computing Machinery. https://doi.org/10.1145/3501714.3501752.
Pfister, Niklas, Stefan Bauer, and Jonas Peters. 2019. “Learning Stable and Predictive Structures in Kinetic Systems.”Proceedings of the National Academy of Sciences 116 (51): 25405–11. https://doi.org/10.1073/pnas.1905688116.
Runge, Jakob, Sebastian Bathiany, Erik Bollt, et al. 2019. “Detecting and Quantifying Causal Associations in Large Nonlinear Time Series Datasets.”Science Advances 5 (11): eaau4996. https://doi.org/10.1126/sciadv.aau4996.
Shi, Jifan, Yang Li, Juan Zhao, et al. 2026. “Deciphering Interventional Dynamical Causality from Non-Intervention Complex Systems.”The Innovation 7 (10): 101358. https://doi.org/10.1016/j.xinn.2026.101358.
Spirtes, Peter, Clark Glymour, and Richard Scheines. 2000. Causation, Prediction, and Search. 2nd ed. MIT Press.
Sugihara, George, Robert May, Hao Ye, et al. 2012. “Detecting Causality in Complex Ecosystems.”Science 338 (6106): 496–500. https://doi.org/10.1126/science.1227079.
Sun, Jie, Dane Taylor, and Erik M. Bollt. 2015. “Causal Entropy Approach to Detect Causal Time Series.”Entropy 17 (11): 7676–702. https://doi.org/10.3390/e17117676.
Zhang, Jiji. 2008. “On the Completeness of Orientation Rules for Causal Discovery in the Presence of Latent Confounders and Selection Bias.”Artificial Intelligence 172 (16): 1873–96. https://doi.org/10.1016/j.artint.2008.08.001.
