Causal discovery: learning the graph from data

Before this lesson: lesson 03 (d-separation, v-structures, the Markov property and faithfulness) and lesson 04 (the do-operator). From probability: conditional independence (probability 04), partial correlation in the multivariate Gaussian (probability 05) and hypothesis tests (probability 11). Time: about 90 minutes, plus the notebook.

---

Why this matters #

• Science. Biologists learn gene regulatory and protein signalling networks from measurements; the classic study of Sachs et al. (2005) recovered most known links of a signalling network from single-cell data with interventions. Climate and Earth scientists use the same ideas on time series (for example the PCMCI method of Runge et al. 2019).
• Engineering and operations. Root-cause analysis in large software systems learns which service's latency drives which other one; discovery on monitoring data suggests the graph that an engineer then checks.
• Machine learning. The adjustment rules of lessons 04–08 need a graph. When nobody can draw it with confidence, discovery gives candidate graphs, and the equivalence class shows which conclusions do not depend on the unknown directions. Discovery ideas also appear in feature selection (Markov blankets) and in invariance methods for distribution shift (lesson 11).
• A lesson in humility. Discovery shows exactly what data alone can and cannot say about causes. This protects you from strong claims based on a learned graph.

1. Predict first #

2. The picture #

2.1 What the data can see: equivalence classes #

Recall the three junctions of lesson 03. In a chain $A \to B \to C$ and in a fork $A \leftarrow B \to C$, the ends $A$ and $C$ are dependent, and independent given $B$. In a collider $A \to B \leftarrow C$ it is the other way round: $A$ and $C$ are independent, and dependent given $B$.

So the chain, the reversed chain and the fork imply the same independence facts. A dataset that fits one fits the other two equally well. The collider is different: it is the only DAG on this skeleton with its independence facts.

Two DAGs are Markov equivalent if they imply the same d-separations, so the same conditional independences. The set of all DAGs that are equivalent to a given one is its Markov equivalence class. The CPDAG (completed partially directed acyclic graph, also called the essential graph) draws a class as one graph: an edge is an arrow if all DAGs in the class orient it the same way, and an undirected line if they disagree.

Interactive figure: loading…

The two arrows of a collider whose ends are not adjacent form a v-structure (also called an immorality: two parents that are "not married"). The figure suggests the main theorem of this lesson, proved in section 3.2: two DAGs are Markov equivalent exactly when they have the same skeleton and the same v-structures.

2.2 The PC algorithm in pictures #

If equivalent graphs cannot be separated, the best we can hope for is the class, drawn as a CPDAG. The PC algorithm (named after its authors Peter Spirtes and Clark Glymour) finds it in three phases:

1. Remove edges. Start with every pair of nodes joined by an undirected edge. For each edge $i - j$, look for a set $S$ of neighbors with $i \perp\!\!\!\perp j \mid S$. If you find one, remove the edge and store $S$ as the separating set $\mathrm{sep}(i, j)$. Try sets of size 0, then 1, then 2, and so on.
2. Orient v-structures. For every triple $i - k - j$ where $i$ and $j$ are not adjacent: if $k$ is not in $\mathrm{sep}(i, j)$, orient $i \to k \leftarrow j$.
3. Propagate. Orient more edges with Meek's rules: an orientation is forced if the other direction would create a new v-structure or a cycle.

Interactive figure: loading…

In the sprinkler graph ($A$ season, $B$ sprinkler, $C$ rain, $D$ wet street, $E$ slippery street, as in lesson 03), PC with oracle tests ends with $A - B$, $A - C$, $B \to D \leftarrow C$ and $D \to E$. Section 4 does every step by hand.

2.3 Breaking the symmetry: the shape of the noise #

Independences cannot orient $A - B$ in the sprinkler graph. Sometimes the shape of the distribution can. Take two variables with $y = x + e$, where the cause $x$ and the noise $e$ are independent. Regress $y$ on $x$: the residual is $e$, which is independent of $x$ by construction. Regress $x$ on $y$: least squares makes the residual uncorrelated with $y$, but if the noise is not Gaussian, the residual is still dependent on $y$. The wrong direction leaves a visible pattern.

Interactive figure: loading…

3. Definitions and main results #

3.1 The question and its four assumptions #

The data are samples from $P$. We want the DAG $G$ of the structural causal model that generated them. Discovery needs a bridge between the graph and the distribution:

1. Causal Markov condition. Every d-separation in $G$ is a conditional independence in $P$: $i \perp_G j \mid S \Rightarrow i \perp\!\!\!\perp j \mid S$. It holds for every SCM with independent noises (lesson 03, section 3.5).
2. Faithfulness. The converse: every conditional independence in $P$ comes from a d-separation, $i \perp\!\!\!\perp j \mid S \Rightarrow i \perp_G j \mid S$. It fails when effects cancel exactly (lesson 03: a direct effect $a$ and an indirect effect $bc$ with $a + bc = 0$).
3. Causal sufficiency. No unmeasured common cause of two measured variables, and no selection of units by a common effect.
4. Acyclicity. The graph has no directed cycles (no feedback loops in the measured time scale).

We also need i.i.d. samples and, on finite data, correct test decisions.

Under assumptions 1–2, the independences of $P$ are exactly the d-separations of $G$. So $P$ tells us the d-separation list of $G$, and nothing more about $G$. Which assumptions can data check? Almost none directly. Faithfulness fails only for special parameter values (a set of measure zero for linear Gaussian and for discrete models), but near-cancellation makes tests weak on finite samples. Causal sufficiency is untestable in general; the FCI algorithm (section 3.3) drops it. Acyclicity is a modeling choice. Some violations leave traces: hidden confounders or cycles can produce independence patterns that no DAG explains, and PC may then report conflicting orientations.

3.2 Markov equivalence #

Definition. DAGs $G_1$ and $G_2$ on the same nodes are Markov equivalent if $i \perp_{G_1} j \mid S \iff i \perp_{G_2} j \mid S$ for all $i, j, S$.

Theorem (Verma and Pearl 1990). Two DAGs are Markov equivalent if and only if they have the same skeleton and the same v-structures.

Two facts about a single DAG explain the theorem. They are also the facts PC uses.

• Fact 1 (adjacency). $i$ and $j$ are adjacent in $G$ if and only if no set $S$ d-separates them. If they are not adjacent, the parents of one of them separate them: $\mathrm{pa}(j)$ if $i$ is not a descendant of $j$, otherwise $\mathrm{pa}(i)$.
• Fact 2 (colliders). Let $i - k - j$ with $i, j$ not adjacent. If $k$ is a collider ($i \to k \leftarrow j$), then no separating set of $i$ and $j$ contains $k$. If $k$ is not a collider, every separating set contains $k$.

Why the two facts hold, and the proof of the theorem

1. If $i \to j$, the one-edge path $i \to j$ has no middle node, so no set blocks it. So adjacent nodes are never d-separated.
2. If $i$ and $j$ are not adjacent and $i$ is not a descendant of $j$, the local Markov property says that $j$ is d-separated from its non-descendants by its parents. So $S = \mathrm{pa}(j)$ separates. (One of the two must be a non-descendant of the other, because $G$ is acyclic.)
3. Fact 2, collider case: the path $i \to k \leftarrow j$ is open whenever $k$ or one of its descendants is in $S$. So a separating set contains neither $k$ nor its descendants.
4. Fact 2, non-collider case: the path $i - k - j$ through a chain or fork is open unless $k \in S$. So every separating set contains $k$.
5. "Only if" of the theorem: by Fact 1 the d-separation list fixes the skeleton. By Fact 2 it fixes, for every unshielded triple, whether the middle node is a collider. So equivalent DAGs share the skeleton and the v-structures.
6. "If" is the harder direction: same skeleton and v-structures imply the same d-separation list. The proof shows that a path is open given $S$ in $G_1$ exactly when a corresponding path is open in $G_2$; see Verma and Pearl (1990) or Pearl (2009), theorem 1.2.8.

Counting. On 3 labeled nodes there are 25 DAGs and 11 equivalence classes:

Skeleton	DAGs	Classes	Class sizes
no edge	1	1	1
one edge (3 choices)	$3 \cdot 2 = 6$	3	2 each ($A \to B$ and $A \leftarrow B$)
two edges (3 choices of middle node)	$3 \cdot 4 = 12$	6	3 (chains and fork) and 1 (collider)
complete	6	1	6 (no v-structure is possible)
total	25	11

The complete graph has 6 DAGs, one per ordering of the 3 nodes: with all three edges present, the data show no independence at all. On 4 nodes there are 543 DAGs in 185 classes; on 10 nodes about $4.2 \cdot 10^{18}$ DAGs. Exhaustive search is impossible beyond a handful of nodes.

3.3 The PC algorithm #

Theorem (Spirtes, Glymour and Scheines). Under causal Markov, faithfulness, causal sufficiency and acyclicity, and with correct independence decisions, the PC algorithm returns the CPDAG of the true DAG.

The three phases, with the reason for each:

PC(variables 1..d, test):
  C := complete undirected graph;  level := 0
  repeat                                      # phase 1: skeleton
    for each edge i - j in C:
      for each S ⊆ adj_C(i) \ {j} with |S| = level:      # also the same with i and j swapped
        if test says i ⊥⊥ j | S:  remove i - j;  sep(i, j) := S;  break
    level := level + 1
  until no node has more than `level` other neighbours
  for each triple i - k - j with i, j not adjacent:       # phase 2: v-structures
    if k ∉ sep(i, j):  orient i → k ← j
  apply Meek's rules R1-R3 until nothing changes         # phase 3

• Phase 1 relies on Fact 1. A non-adjacent pair is separated by the parents of one of them, and the parents are among the neighbors, so searching subsets of the current neighbors is enough. It starts with small sets because they are cheap and their tests are more reliable.
• Phase 2 relies on Fact 2. One separating set is enough: if $k$ is not in it, $k$ is a collider.
• Phase 3 orients the edges that every DAG of the class shares:

Rule	If you see	Then orient	Otherwise you would create
R1	$a \to b - c$, with $a, c$ not adjacent	$b \to c$	a new v-structure $a \to b \leftarrow c$
R2	$a \to c \to b$ and $a - b$	$a \to b$	a cycle $a \to c \to b \to a$
R3	$a - c \to b$, $a - e \to b$, $a - b$, with $c, e$ not adjacent	$a \to b$	a cycle or a new v-structure (check both cases of $b \to a$)

Meek (1995) proved that R1–R3 orient exactly the edges that are directed in the CPDAG. (A fourth rule is needed only when you add background knowledge, such as "$B$ cannot cause $A$".)

The version used here. In the original PC, the result of phase 1 can depend on the order in which pairs are tested, because removing one edge changes the neighbor sets of later tests. PC-stable (Colombo and Maathuis 2014) freezes the neighbor sets at the start of each level, so the order no longer matters. The widget and the notebook use PC-stable.

Cost. The number of tests grows like $d^{q}$, where $q$ is the largest number of neighbors of a node. PC is fast on sparse graphs with hundreds of variables and slow on dense ones.

Independence tests on data. For (approximately) Gaussian data, $i \perp\!\!\!\perp j \mid S$ is equivalent to a zero partial correlation $\rho_{ij \cdot S}$. With $P$ the inverse of the correlation matrix of the variables $(i, j, S)$,

Under independence, $z$ is approximately $\mathcal N(0, 1)$ (Fisher's z test). PC accepts independence when the p-value is larger than the level $\alpha$. For nonlinear or discrete data, use a test of conditional mutual information or a kernel test (link: information theory 07).

Hidden common causes: FCI. If causal sufficiency fails, PC can return an edge $X - Y$ for two variables that only share a hidden cause. The FCI algorithm (Fast Causal Inference, Spirtes et al.) allows hidden common causes and selection. Its output, a partial ancestral graph, marks each edge end as an arrowhead, a tail or a circle ("cannot decide"); $X \leftrightarrow Y$ means "a hidden common cause, and no causal effect either way". It needs more tests and usually orients fewer edges.

3.4 Score-based methods and continuous optimization #

Instead of tests, give every graph a score and search for the best one. For linear-Gaussian or discrete models a standard score is the BIC:

where the first term is the maximized log-likelihood of each variable given its parents and $k_G$ the number of parameters. The penalty prefers sparse graphs. Markov-equivalent DAGs get the same BIC (score equivalence), so scores also cannot go beyond the class.

Greedy equivalence search (GES, Chickering 2002) moves between equivalence classes: a forward phase adds edges while the score improves, a backward phase removes them. With a consistent, score-equivalent score, GES finds the true class as $n \to \infty$ (under faithfulness), although it looks at only a tiny part of the search space.

NOTEARS (Zheng et al. 2018) turns the combinatorial search into a smooth problem. For a linear SCM $X = XW + N$ with a weighted adjacency matrix $W \in \R^{d \times d}$:

where $W \circ W$ squares each entry and $e^{A} = \sum_{k \ge 0} A^k / k!$ is the matrix exponential.

Theorem. $h(W) \ge 0$, and $h(W) = 0$ if and only if the graph of $W$ (an edge $i \to j$ when $W_{ij} \neq 0$) is acyclic.

Proof

1. Let $A = W \circ W$. All entries of $A$ are $\ge 0$, and $A_{ij} > 0$ exactly when $i \to j$ is an edge.
2. $(A^k)_{ii}$ is a sum over all closed walks $i \to \dots \to i$ of length $k$ of the products of the $A$ entries along the walk. Each product is $\ge 0$, and it is $> 0$ when all edges of the walk exist. So $\operatorname{tr}(A^k) > 0$ if the graph has a closed walk of length $k$, and $= 0$ otherwise.
3. $\operatorname{tr}(e^A) = \operatorname{tr}(I) + \sum_{k \ge 1} \operatorname{tr}(A^k)/k! = d + \sum_{k \ge 1} \operatorname{tr}(A^k)/k!$. So $h(W) = \sum_{k \ge 1} \operatorname{tr}(A^k)/k! \ge 0$.
4. A graph has a directed cycle if and only if it has a closed walk of some length $k \ge 1$ (a self-loop is the case $k = 1$). So $h(W) = 0$ exactly for acyclic graphs. $\blacksquare$

Small example. For two nodes with $W = \begin{pmatrix} 0 & a \\ b & 0 \end{pmatrix}$ (edges $1 \to 2$ with weight $a$, $2 \to 1$ with weight $b$), $A^2 = a^2 b^2 I$ and the odd powers have zero trace, so $h(W) = 2\cosh(ab) - 2$. It is zero exactly when $a = 0$ or $b = 0$. For $a = 1$, $b = 0.5$: $h = 2\cosh(0.5) - 2 \approx 0.255$. The gradient $\nabla h(W) = (e^{W \circ W})^\top \circ 2W$ lets a solver (an augmented Lagrangian method in the paper) push $h$ to zero. The problem is not convex, so the solver finds a local solution, and small weights are set to zero at the end.

The varsortability warning. NOTEARS and its successors did very well on simulated benchmarks. Reisach, Seiler and Weichwald (2021) found a reason that has nothing to do with causal reasoning. In the usual simulations (weights with $|w| \in [0.5, 2]$, equal noise variances), the marginal variance grows along the causal order: each effect adds its own noise to the variance it receives from its causes. A trivial baseline, "sort the variables by variance, then regress each one on the earlier ones", is as good as NOTEARS there. After standardizing each variable to variance 1, these methods get much worse, although standardizing changes no causal fact.

Interactive figure: loading…

The practical lesson: a discovery method should give the same graph when you change the units of a variable (metres or kilometres). Check this by running it on standardized data too.

3.5 Beyond equivalence: functional assumptions #

Markov equivalence is a limit for methods that use only conditional independences. A method that assumes more about the functions or the noise can sometimes orient every edge.

LiNGAM (Linear Non-Gaussian Acyclic Model; Shimizu et al. 2006): $X = BX + e$ with $B$ the weights of a DAG and $e_1, \dots, e_d$ independent and non-Gaussian. Then $B$, and so the full DAG, is identifiable from $P$. The original algorithm uses independent component analysis (ICA): $X = (I - B)^{-1}e$ is a mixture of independent non-Gaussian sources, and ICA can recover the mixing matrix.

Two variables, why it works. Let $y = x + e$ with $x \perp\!\!\!\perp e$, both with variance 1. The backward regression of $x$ on $y$ has slope $b = \operatorname{Cov}(x, y)/\operatorname{Var}(y) = 1/2$, so its residual is

Both $r$ and $y$ are linear combinations of the independent variables $x$ and $e$, with all four coefficients nonzero. The Darmois–Skitovich theorem says: if two such linear combinations are independent, then every variable that appears in both is Gaussian. So if $x$ or $e$ is not Gaussian, $r$ and $y$ are dependent, although they are uncorrelated. In the causal direction the residual is $e$, which is independent of $x$. The asymmetry identifies the direction. With uniform noise you can see the dependence: when $y$ is near its largest value $2\sqrt3$, both $x$ and $e$ must be near $\sqrt3$, so $r = (x - e)/2$ must be near 0. This is the diamond in the figure of section 2.3.

Additive noise models (Hoyer et al. 2009) extend the idea to nonlinear functions: $y = f(x) + e$ with $e \perp\!\!\!\perp x$. For generic nonlinear $f$, no model $x = g(y) + \tilde e$ with $\tilde e \perp\!\!\!\perp y$ exists. The test is the same: fit both directions with a regression method and keep the direction whose residual is independent of the input. The linear-Gaussian case is one of the few exceptions where both directions fit.

Gaussian noise with equal variances is identifiable too

Peters and Bühlmann (2014) proved that a linear-Gaussian SCM is identifiable if all noise variances are equal (the common value $\sigma^2$ may be unknown). One way to see it: take the variables in some order and compute the conditional variance of each one given the earlier ones. In a causal order each of these equals $\sigma^2$. In every order their product is $\det \Sigma$, so by the inequality of arithmetic and geometric means their sum is at least $d\,(\det\Sigma)^{1/d} = d\sigma^2$, with equality only when they are all equal. The theorem says that a non-causal order cannot make them all equal. So the orders with the smallest sum of conditional variances are exactly the causal orders. This assumption is strong and not invariant to rescaling. It is closely related to the varsortability warning: equal noise variances in the raw units are exactly what makes simulated variances grow along the order.

3.6 Interventions shrink the class #

An experiment changes the picture. A perfect intervention $\mathrm{do}(B = b)$ cuts all arrows into $B$ (lesson 04). Compare the observational and the interventional data:

• If $A \to B$, setting $B$ by hand does not change the distribution of $A$.
• If $B \to A$, setting $B$ changes the distribution of $A$.

So one experiment on $B$ orients every edge at $B$. Then Meek's rules may orient more. Two DAGs that also agree on the edges at every intervened node form an interventional equivalence class (Hauser and Bühlmann 2012). Eberhardt, Glymour and Scheines (2005) showed that $d - 1$ experiments on single variables are always enough to orient every edge among $d$ variables, and in the worst case necessary.

Interactive figure: loading…

3.7 Honest status of causal discovery #

• Hidden variables, measurement error and selection break causal sufficiency and the meaning of the tests. Most real datasets have all three.
• Tests are weak for conditioning sets that are large compared with $n$, and for nonlinear dependence; one wrong test can change many edges.
• Different methods disagree on real data, and benchmark results on simulated data can be misleading (varsortability).
• Good use: generate hypotheses, combine the output with domain knowledge (time order, known absences of effects), report the equivalence class rather than one DAG, check the stability of the edges under resampling (bootstrap), and confirm important edges with experiments.

Answers to the predictions. (1) No. The chain and the reversed chain have the same skeleton and no v-structure, so they imply the same independences; no amount of observational data separates them without extra assumptions. (2) The collider $A \to B \leftarrow C$: it is the only graph in which $A$ and $C$ are independent but dependent given $B$. (3) Yes. With uniform noise the residual of the anti-causal regression is uncorrelated with $y$ but not independent of it (section 3.5).

4. Worked example #

Run PC with oracle tests on the sprinkler graph: $A \to B$, $A \to C$, $B \to D$, $C \to D$, $D \to E$ ($A$ season, $B$ sprinkler, $C$ rain, $D$ wet, $E$ slippery). The oracle answers each test with d-separation in the true graph.

Phase 1, level 0 ($S = \varnothing$). Test all $\binom52 = 10$ pairs for marginal independence. Every pair is joined by an open path (for example $B \leftarrow A \to C$ for $B, C$ and $A \to B \to D \to E$ for $A, E$). No edge is removed.

Phase 1, level 1 ($|S| = 1$). All neighbor sets are still complete. The tests that succeed:

Pair	Separating set	Why the paths are blocked
$A, E$	$\{D\}$	every path from $A$ to $E$ ends with $D \to E$, and $D$ is a non-collider in $S$
$B, C$	$\{A\}$	$B \leftarrow A \to C$ is blocked by $A$; $B \to D \leftarrow C$ is blocked at the collider $D$
$B, E$	$\{D\}$	every path ends with $D \to E$
$C, E$	$\{D\}$	the same

$A, D$ survives level 1: given $B$ the path $A \to C \to D$ is open, and given $C$ the path $A \to B \to D$ is open. Given $E$ both are open. After level 1 the graph is $A - B$, $A - C$, $A - D$, $B - D$, $C - D$, $D - E$.

Phase 1, level 2. Now $\mathrm{adj}(A) \setminus \{D\} = \{B, C\}$. The test $A \perp\!\!\!\perp D \mid \{B, C\}$ succeeds: both directed paths are blocked, and there is no other path. Remove $A - D$ with $\mathrm{sep}(A, D) = \{B, C\}$. The remaining five edges are true edges; no set separates them. Level 3 is not possible (no node has 4 other neighbors left), so phase 1 ends. The skeleton is now exactly the true skeleton.

Phase 2. Check every unshielded triple (two ends not adjacent):

Triple	$\mathrm{sep}$ of the ends	Middle node in it?	Result
$B - A - C$	$\{A\}$	yes	no v-structure
$B - D - C$	$\{A\}$	no	$B \to D \leftarrow C$
$A - B - D$	$\{B, C\}$	yes	nothing
$A - C - D$	$\{B, C\}$	yes	nothing
$B - D - E$, $C - D - E$	$\{D\}$	yes	nothing

Phase 3. R1 on $B \to D - E$ ($B$ and $E$ not adjacent) gives $D \to E$. No arrow points into $A$, $B$ or $C$, so no rule applies to $A - B$ and $A - C$.

Result. The CPDAG is $A - B$, $A - C$, $B \to D \leftarrow C$, $D \to E$. Its class has 3 DAGs:

1. $A \to B$, $A \to C$ (the truth);
2. $B \to A \to C$;
3. $C \to A \to B$.

The fourth orientation, $B \to A \leftarrow C$, is not in the class: it adds a v-structure at $A$ ($B$ and $C$ are not adjacent). The data cannot say whether the season drives the sprinkler or the other way round. Background knowledge does: nothing in the graph can cause the season, so arrows into $A$ are impossible, and only DAG 1 remains.

PC needed only 5 successful tests for the removals, but it ran 78 tests in total (the widget counts them). On data, each of the 78 can be wrong.

5. In code #

NumPy and SciPy count the equivalence classes, run the key Fisher-z tests of the worked example on simulated linear-Gaussian data, and test the two directions of $y = x + e$.

import itertools
import numpy as np
from scipy import stats

# (a) count DAGs and Markov equivalence classes on 3 labelled nodes
pairs = list(itertools.combinations(range(3), 2))
def acyclic(E):                      # E: set of directed edges (i, j)
    return not any((a, b) in E and (b, c) in E and (c, a) in E for a, b, c in itertools.permutations(range(3)))
def key(E):                          # skeleton + v-structures
    skel = frozenset(frozenset(e) for e in E)
    vs = frozenset((frozenset((i, j)), k) for (i, k), (j, k2) in itertools.permutations(E, 2)
                   if k == k2 and i != j and frozenset((i, j)) not in skel)
    return skel, vs
dags = []
for s in itertools.product((0, 1, 2), repeat=3):
    E = {(i, j) if v == 1 else (j, i) for (i, j), v in zip(pairs, s) if v}
    if acyclic(E):
        dags.append(E)
print(len(dags), len({key(E) for E in dags}))                       # 25 11

# (b) Fisher-z tests on the sprinkler SCM: A season, B sprinkler, C rain, D wet, E slippery
W = np.zeros((5, 5))
for (i, j), w in {(0, 1): 0.8, (0, 2): 0.7, (1, 3): 0.9, (2, 3): 0.8, (3, 4): 0.7}.items():
    W[i, j] = w
rng = np.random.default_rng(0)
n = 1000
X = np.zeros((n, 5))
for j in range(5):                                                   # nodes are in causal order
    X[:, j] = X @ W[:, j] + rng.normal(size=n)
C = np.corrcoef(X.T)

def fisher_z(i, j, S):
    idx = [i, j, *S]
    P = np.linalg.inv(C[np.ix_(idx, idx)])
    r = -P[0, 1] / np.sqrt(P[0, 0] * P[1, 1])                       # partial correlation
    z = np.arctanh(r) * np.sqrt(n - len(S) - 3)
    return r, 2 * stats.norm.sf(abs(z))

for i, j, S in [(1, 2, []), (1, 2, [0]), (0, 4, [3]), (0, 3, [1]), (0, 3, [1, 2])]:
    r, p = fisher_z(i, j, S)
    print("ABCDE"[i], "ABCDE"[j], ["ABCDE"[s] for s in S], round(r, 3), round(p, 3))
# B C [] 0.35 0.0
# B C ['A'] 0.008 0.809
# A E ['D'] -0.055 0.08
# A D ['B'] 0.289 0.0
# A D ['B', 'C'] 0.001 0.968

# (c) the direction of y = x + e with uniform noise: residual independent of the regressor?
def hsic(a, b):
    K = np.exp(-(a[:, None] - a[None]) ** 2 / np.median((a[:, None] - a[None]) ** 2))
    L = np.exp(-(b[:, None] - b[None]) ** 2 / np.median((b[:, None] - b[None]) ** 2))
    H = np.eye(len(a)) - 1 / len(a)
    return np.sum((H @ K @ H) * L) / len(a) ** 2
x = rng.uniform(-np.sqrt(3), np.sqrt(3), 400)
y = x + rng.uniform(-np.sqrt(3), np.sqrt(3), 400)
res_y = y - np.polyval(np.polyfit(x, y, 1), x)                       # fit y from x
res_x = x - np.polyval(np.polyfit(y, x, 1), y)                       # fit x from y
print(round(np.corrcoef(x, res_y)[0, 1], 3), round(np.corrcoef(y, res_x)[0, 1], 3))   # 0.0 0.0
print(round(1000 * hsic(x, res_y), 2), round(1000 * hsic(y, res_x), 2))              # 0.62 6.54

The tests agree with the worked example: $B$ and $C$ are dependent but independent given $A$, and $A$ and $D$ are independent only given both $B$ and $C$. The test of $A, E$ given $D$ has $p = 0.08$: with $\alpha = 0.1$ it would count as dependent. (PC still recovers the correct CPDAG on this sample for every $\alpha$ from 0.01 to 0.2: when this test keeps $A - E$, another test with a different separating set removes it.) In (c) both residuals are uncorrelated with their regressor, but the HSIC statistic is ten times larger in the anti-causal direction. The notebook builds the full PC algorithm and a permutation test.

6. Where it shows up #

• Molecular biology. Sachs et al. (2005) learned a protein signalling network in human immune cells from flow-cytometry data with interventions; it is still a standard test case. Gene regulatory network inference uses constraint-based and score-based methods on expression data, plus knock-out experiments as interventions.
• Climate and Earth science. PCMCI (Runge et al. 2019) adapts PC to time series, where time order fixes most directions, to find drivers of climate patterns such as monsoons and El Niño teleconnections.
• Neuroscience. Effective connectivity between brain regions is estimated from fMRI and electrophysiology data with discovery methods, with all the caveats of hidden variables and measurement error.
• Industry. Root-cause analysis in microservice systems and manufacturing lines learns a graph of metrics, then follows the arrows back from an anomaly. Feature selection uses the Markov blanket (parents, children and the children's other parents) of the target, which PC-like algorithms find.
• Causal ML. Discovery gives candidate graphs for the adjustment of lessons 04–08; the equivalence class shows whether an adjustment set is valid in every graph of the class.

7. Common mistakes #

• Reading every arrow of a learned graph as causal. Undirected edges of a CPDAG can point either way, and an arrow is only as good as the tests and the assumptions behind it. Report the class and say which assumptions you made.
• Forgetting hidden confounders. PC assumes causal sufficiency. With a hidden common cause it draws an edge between two variables that have no effect on each other. Use FCI when hidden confounding is plausible, and check whether edges survive.
• Treating "not rejected" as "independent" without power. A small $n$ removes true edges. Check the stability of each edge under bootstrap resampling and with several values of $\alpha$.
• Using the wrong test. Fisher's z test only sees linear dependence. For nonlinear or discrete data, use a kernel or conditional mutual information test, or the result can miss real edges.
• Trusting results that change with the units. If standardizing the data changes the graph, the method used the scale (varsortability), not the causal structure.
• Ignoring time order and domain knowledge. A cause cannot come after its effect. Add such constraints as background knowledge; they orient edges that the data cannot.
• Selecting units by an outcome before discovery. Selection on a common effect creates dependences (Berkson's paradox, lesson 03) that discovery then reads as edges.

8. Check yourself #

1. How many DAGs are in the Markov equivalence class of the chain $A \to B \to C \to D$? Draw its CPDAG.
   Answer

   4. The CPDAG is $A - B - C - D$ with all edges undirected. A DAG in the class must have no v-structure, so at most one arrow can point into each node from the path; such an orientation is fixed by choosing the one node where all arrows start (the source): $A$, $B$, $C$ or $D$.
2. Why does phase 1 of PC only test conditioning sets made of current neighbors of $i$ or of $j$?
   Answer

   If $i$ and $j$ are not adjacent in the true DAG, the parents of one of them separate them (Fact 1). Parents are neighbors in the true graph, and the current graph of PC always contains the true skeleton (with correct tests), so some subset of the current neighbors is enough.
3. Compute $h(W) = \operatorname{tr}(e^{W \circ W}) - d$ for $W = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}$ and for $W = \begin{pmatrix} 0 & 1 \\ 0.5 & 0 \end{pmatrix}$.
   Answer

   First: $A = W \circ W$ is nilpotent ($A^2 = 0$), so $e^A = I + A$, trace 2, $h = 0$: a DAG. Second: $A^2 = 0.25\,I$, so $\operatorname{tr}(e^A) = 2\cosh(0.5) \approx 2.255$ and $h \approx 0.255 > 0$: the graph has the cycle $1 \to 2 \to 1$.
4. In the sprinkler class (section 4), which DAG is left if you know that nothing in the graph causes the season? Without this knowledge, how many DAGs are left after one experiment on $B$?
   Answer

   With the background knowledge only the true DAG is left, because the other two members have an arrow into $A$. The experiment on $B$ orients $A \to B$ (and $B \to D$, already oriented). No Meek rule then applies to $A - C$, because no arrow points into $A$. So $A - C$ stays undirected and 2 DAGs are left: $A \to C$ and $C \to A$.
5. Two variables $X$ and $Y$ have a hidden common cause and no effect on each other. What does PC return, and what does FCI return?
   Answer

   PC returns the edge $X - Y$: they are dependent and there is no other variable to condition on. FCI returns $X$ o–o $Y$ (with only two variables it cannot decide); with suitable other variables it can return $X \leftrightarrow Y$, which says "hidden common cause, no direct effect".
6. A colleague runs NOTEARS on raw sensor data and gets a clean chain. On standardized data the chain is reversed. What do you conclude?
   Answer

   The result depends on the scale, so it is probably driven by the variances (varsortability), not by causal structure. Units are a convention and cannot carry causal information unless a strong assumption (such as equal noise variances in these units) is justified. Report the result as unreliable and try scale-invariant methods (PC with partial correlations, GES with BIC) or experiments.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What does faithfulness assume?

   Show answer

2. When are two DAGs Markov equivalent?

   Show answer

3. What is a CPDAG?

   Show answer

4. How many DAGs and Markov equivalence classes are there on 3 labeled nodes?

   Show answer

5. What are the three phases of the PC algorithm?

   Show answer

6. In PC, when is $i - k - j$ (with $i, j$ not adjacent) oriented as $i \to k \leftarrow j$?

   Show answer

7. What is the NOTEARS acyclicity function, and why is it zero exactly for DAGs?

   Show answer

8. Why can LiNGAM orient an edge that PC cannot?

   Show answer

9. What is the varsortability warning?

   Show answer

10. Why does an intervention on $B$ help discovery?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Can you learn causal structure from observational data?" — Partly. Under causal Markov, faithfulness, no hidden confounders and acyclicity, you can learn the Markov equivalence class (skeleton plus v-structures, then Meek's rules), not a single DAG. Orienting the rest needs extra assumptions (non-Gaussian or nonlinear additive noise), background knowledge or experiments.
• "How does the PC algorithm work, and what can go wrong?" — Three phases (skeleton by independence tests, v-structures from separating sets, Meek's rules). On finite data, a missed dependence deletes a true edge, a wrong separating set creates a wrong v-structure, and errors propagate. Hidden confounders produce spurious edges; FCI handles them.
• "Would you trust NOTEARS on our data?" — Only with care: the problem is non-convex, the least squares score depends on the scale, and benchmark wins can come from varsortability. Run it on standardized data, compare with PC or GES, check stability under bootstrap, and validate key edges with experiments or domain knowledge.
• "A learned graph says feature X causes churn. What next?" — Treat it as a hypothesis. Check the assumptions (hidden confounders, selection), the stability of the edge, and whether its direction is decided by the data or only by the method; then test it with an A/B experiment.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. Write fisher_z_test(C, n, i, j, S): partial correlation from the correlation matrix and the p-value of Fisher's z test.
2. Write pc_skeleton and orient_v_structures (Meek's rules are given) and measure the structural Hamming distance to the true CPDAG of a 5-node linear-Gaussian SCM for several $n$ and $\alpha$.
3. Enumerate the Markov equivalence class of a small DAG by brute force and check the counts 25 and 11 on 3 nodes.
4. Test the direction of $y = x + e$ with non-Gaussian noise: regress both ways and compare a simple HSIC permutation test of residual and regressor.
5. Varsortability: show that "sort by variance" finds the causal order on standard simulated data and fails after standardizing.

11. Watch and read #

10 - Causal Discovery from Observational DataBrady Neal - Causal Inference · 51 min

The full lecture: assumptions, Markov equivalence, PC, and the semi-parametric (non-Gaussian, additive noise) setting.

10.4 - The PC Algorithm for Causal DiscoveryBrady Neal - Causal Inference · 14 min

The three phases of PC on a small example, as in section 4.

10.6 - Linear Non-Gaussian Setting for Causal DiscoveryBrady Neal - Causal Inference · 7 min

Why non-Gaussian noise identifies the direction of a linear edge.

11 - Causal Discovery from InterventionsBrady Neal - Causal Inference · 51 min

How experiments shrink the equivalence class, and how many are needed.

LiNGAM approach to causal discovery. Talk at the KDD2021 Workshop on Causal Discovery (CD2021)SHIMIZU Shohei · 56 min

LiNGAM and its extensions, presented by its first author.

Causality, part 1 - Bernhard Schölkopf - MLSS 2020, Tübingenvirtual mlss2020 · 1 h 36 min

The machine-learning view of causal models and why independence of mechanisms helps discovery.

• Peters, Janzing and Schölkopf, Elements of Causal Inference (MIT Press, 2017, open access): chapter 4 (two-variable identifiability: LiNGAM and additive noise) and chapter 7 (learning multivariate models: PC, GES, interventions).
• Brady Neal, Introduction to Causal Inference (free PDF), chapters 11–12 (causal discovery from observational and interventional data).
• Glymour, Zhang and Spirtes (2019), Review of causal discovery methods based on graphical models, Frontiers in Genetics (open access).
• Colombo and Maathuis (2014), Order-independent constraint-based causal structure learning (PC-stable).
• Chickering (2002), Optimal structure identification with greedy search (GES).
• Zheng, Aragam, Ravikumar and Xing (2018), DAGs with NO TEARS: continuous optimization for structure learning.
• Reisach, Seiler and Weichwald (2021), Beware of the simulated DAG! Causal discovery benchmarks may be easy to game.
• Shimizu, Hoyer, Hyvärinen and Kerminen (2006), A linear non-Gaussian acyclic model for causal discovery (LiNGAM).
• Hoyer, Janzing, Mooij, Peters and Schölkopf (2009), Nonlinear causal discovery with additive noise models.
• Peters and Bühlmann (2014), Identifiability of Gaussian structural equation models with equal error variances.
• Hauser and Bühlmann (2012), Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs.