Causal Inference VIII(b) — Causal Discovery: learning the DAG from data¶
PC, GES, and LiNGAM — what observational data can (and cannot) tell you about causal structure¶
The DAG notebook assumed the causal graph was known and read identification off it. But where does the graph come from? Domain knowledge, usually — yet a natural question is whether the structure itself can be learned from data. That is causal discovery, and its answer is subtle and important: from purely observational data you can recover the graph's skeleton and some edge directions, but generally only up to a Markov equivalence class — a set of DAGs that imply the same conditional independencies and are therefore statistically indistinguishable. Extra assumptions (or non-Gaussianity) are needed to orient the rest.
This notebook builds the two main families and the sharp result that separates them:
- Constraint-based (PC algorithm) — test conditional independencies to find the skeleton, then orient v-structures (unshielded colliders) and propagate with logical rules. It returns a CPDAG (completed partially directed acyclic graph): some edges directed, some left undirected because the data cannot distinguish their orientation.
- Score-based (GES) — greedily search over graphs to maximize a fit score (BIC); it also returns a CPDAG.
- LiNGAM (Linear Non-Gaussian Acyclic Model) — the key insight (Shimizu et al. 2006): if the noise is non-Gaussian, the causal direction becomes identifiable (the regression residuals are independent of the cause only in the true direction), so LiNGAM recovers the entire DAG, orienting edges the CPDAG methods cannot.
We first simulate from a known DAG (the only way to prove an algorithm recovers the truth) and show PC leaving an edge undirected while LiNGAM recovers every direction; then we apply discovery to real single-cell data — Sachs et al.'s (2005) protein-signaling network, the field's benchmark, which has a biologically established ground-truth DAG to score against. Python-lead (causal-learn); R companion uses bnlearn.
1. The PC algorithm — skeleton plus what colliders reveal¶
We generate data from the DAG $X_4\to X_0\to X_2$, $X_1\to X_2$, $X_2\to X_3$ with a linear structural model and Gaussian noise. The PC algorithm first recovers the skeleton (which variables are directly connected) via conditional-independence tests, then orients edges it logically can (colliders and propagation). With Gaussian data it recovers the skeleton and the collider but leaves one edge undirected — the Markov-equivalence limit.
import numpy as np, pandas as pd, matplotlib.pyplot as plt, networkx as nx, warnings
warnings.filterwarnings("ignore")
from causallearn.search.ConstraintBased.PC import pc
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
def sem(n, noise, seed=1):
r=np.random.default_rng(seed); e=noise(r,(n,5))
X4=e[:,4]; X0=0.8*X4+e[:,0]; X1=e[:,1]; X2=0.7*X0+0.6*X1+e[:,2]; X3=0.9*X2+e[:,3]
return np.column_stack([X0,X1,X2,X3,X4]) # true edges: 4->0, 0->2, 1->2, 2->3
lbl={0:"X0",1:"X1",2:"X2",3:"X3",4:"X4"}
pos={0:(1,1),1:(2,0),2:(2,1),3:(3,1),4:(0,1)}
def draw(ax, directed, undirected, title):
for n,(x,y) in pos.items():
ax.scatter([x],[y],s=1400,facecolor="white",edgecolor="k",zorder=3); ax.text(x,y,lbl[n],ha="center",va="center",zorder=4)
for u,v in directed:
x1,y1=pos[u]; x2,y2=pos[v]; dx,dy=x2-x1,y2-y1; L=np.hypot(dx,dy); ux,uy=dx/L,dy/L; r=0.17
ax.annotate("",xy=(x2-ux*r,y2-uy*r),xytext=(x1+ux*r,y1+uy*r),arrowprops=dict(arrowstyle="-|>",color="k",lw=2),zorder=2)
for u,v in undirected:
x1,y1=pos[u]; x2,y2=pos[v]; dx,dy=x2-x1,y2-y1; L=np.hypot(dx,dy); ux,uy=dx/L,dy/L; r=0.17
ax.plot([x1+ux*r,x2-ux*r],[y1+uy*r,y2-uy*r],color=RED,lw=2.5,ls="--",zorder=2)
ax.set_title(title,fontsize=11); ax.axis("off"); ax.set_xlim(-0.6,3.6); ax.set_ylim(-0.5,1.5)
Xg=sem(5000, lambda r,s:r.normal(0,1,s))
cg=pc(Xg, 0.05, show_progress=False); A=cg.G.graph
di=[]; un=[]
for i in range(5):
for j in range(i+1,5):
if A[i,j]==-1 and A[j,i]==1: di.append((i,j))
elif A[i,j]==1 and A[j,i]==-1: di.append((j,i))
elif A[i,j]==-1 and A[j,i]==-1: un.append((i,j))
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
draw(ax[0],[(4,0),(0,2),(1,2),(2,3)],[],"TRUE DAG (the data-generating structure)")
draw(ax[1],di,un,"PC learns the CPDAG (Gaussian data)")
plt.tight_layout(); plt.show()
print(f"PC directed edges: {[f'{lbl[u]}->{lbl[v]}' for u,v in di]}")
print(f"PC undirected (unidentified) edges: {[f'{lbl[u]}-{lbl[v]}' for u,v in un]}")
print("PC recovered the full skeleton, oriented the collider X0->X2<-X1 and propagated X2->X3, but left X0-X4 undirected:")
print("Gaussian observational data identifies only the Markov equivalence class, not every causal direction.")
PC directed edges: ['X0->X2', 'X1->X2', 'X2->X3'] PC undirected (unidentified) edges: ['X0-X4'] PC recovered the full skeleton, oriented the collider X0->X2<-X1 and propagated X2->X3, but left X0-X4 undirected: Gaussian observational data identifies only the Markov equivalence class, not every causal direction.
2. LiNGAM — non-Gaussianity identifies the full DAG¶
The undirected edge is not a failure of the algorithm; it is a fundamental limit of what Gaussian observational data contains. LiNGAM escapes it with one extra ingredient: non-Gaussian noise. The intuition is an asymmetry invisible under Gaussianity — if $Y=\beta X+\varepsilon$ with non-Gaussian $\varepsilon$, the residual is independent of the regressor only in the true causal direction. Re-running on the same structure with uniform noise, LiNGAM recovers every direction, including the one PC left ambiguous.
from causallearn.search.FCMBased import lingam
Xn=sem(5000, lambda r,s:r.uniform(-1.7,1.7,s))
model=lingam.DirectLiNGAM(); model.fit(Xn); B=model.adjacency_matrix_
edges=[(j,i) for i in range(5) for j in range(5) if abs(B[i,j])>0.15]
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
draw(ax[0],[(4,0),(0,2),(1,2),(2,3)],[],"TRUE DAG")
draw(ax[1],edges,[],"LiNGAM recovers the FULL DAG (non-Gaussian data)")
plt.tight_layout(); plt.show()
print(f"LiNGAM causal order: {[lbl[i] for i in model.causal_order_]}")
print(f"LiNGAM edges: {[f'{lbl[u]}->{lbl[v]}' for u,v in edges]}")
correct=set(edges)=={(4,0),(0,2),(1,2),(2,3)}
print(f"Matches the true DAG exactly: {correct}. The X0-X4 edge PC left undirected is now oriented X4->X0 -- correctly.")
print("Non-Gaussianity is the extra structure that breaks Markov-equivalence ties and identifies the whole graph.")
print()
nz=sorted(abs(B[i,j]) for i in range(5) for j in range(5) if i!=j and abs(B[i,j])>0)
print(f"The 0.15 cutoff above looks like a tuning knob and is not one. DirectLiNGAM prunes to exact")
print(f"zeros: the fitted matrix has {len(nz)} non-zero entries, the smallest {min(nz):.2f} and the largest {max(nz):.2f},")
print(f"with nothing in between and nothing below. Any cutoff from 0.05 to 0.40 returns the same graph.")
print(f"Recovered path coefficients {[f'{v:.2f}' for v in nz]} against true 0.6, 0.7, 0.8, 0.9.")
LiNGAM causal order: ['X4', 'X0', 'X1', 'X2', 'X3'] LiNGAM edges: ['X4->X0', 'X0->X2', 'X1->X2', 'X2->X3'] Matches the true DAG exactly: True. The X0-X4 edge PC left undirected is now oriented X4->X0 -- correctly. Non-Gaussianity is the extra structure that breaks Markov-equivalence ties and identifies the whole graph. The 0.15 cutoff above looks like a tuning knob and is not one. DirectLiNGAM prunes to exact zeros: the fitted matrix has 4 non-zero entries, the smallest 0.63 and the largest 0.91, with nothing in between and nothing below. Any cutoff from 0.05 to 0.40 returns the same graph. Recovered path coefficients ['0.63', '0.70', '0.79', '0.91'] against true 0.6, 0.7, 0.8, 0.9.
3. The assumptions — and why discovery is a hypothesis generator, not a design¶
Causal discovery is powerful but rests on strong, largely untestable assumptions:
- Causal sufficiency — no unmeasured common causes. A hidden confounder of two measured variables is mistaken for a direct edge.
- Faithfulness — the only independencies in the data are those implied by the graph (no exact cancellations).
- Reliable conditional-independence tests / correct score — hard in finite samples and high dimensions.
Score-based GES returns the same CPDAG as PC here; and we demonstrate the causal-sufficiency failure explicitly — a hidden common cause of two otherwise-unconnected variables produces a phantom edge.
from causallearn.search.ScoreBased.GES import ges
gr=ges(Xg) # score-based (BIC) on the Gaussian data -> also a CPDAG
Ag=gr["G"].graph; dg=[]; ug=[]
for i in range(5):
for j in range(i+1,5):
if Ag[i,j]==-1 and Ag[j,i]==1: dg.append((i,j))
elif Ag[i,j]==1 and Ag[j,i]==-1: dg.append((j,i))
elif Ag[i,j]==-1 and Ag[j,i]==-1: ug.append((i,j))
print("GES (score-based, Gaussian) also returns a CPDAG:")
print(f" directed: {[f'{lbl[u]}->{lbl[v]}' for u,v in dg]} undirected: {[f'{lbl[u]}-{lbl[v]}' for u,v in ug]}")
r=np.random.default_rng(3); n=5000; H=r.normal(0,1,n) # hidden common cause of A and B
Ahid=H+r.normal(0,1,n); Bhid=H+r.normal(0,1,n) # A and B have NO direct edge
cgh=pc(np.column_stack([Ahid,Bhid]),0.05,show_progress=False)
phantom = cgh.G.graph[0,1]!=0
print(f"\nCausal-sufficiency failure: A and B share a HIDDEN cause H (no true A-B edge).")
print(f" PC (not seeing H) infers a direct edge between A and B: {phantom} -> a phantom edge from unmeasured confounding.")
print("This is why discovery needs causal sufficiency (or FCI), and why a learned edge is a hypothesis to test, not a fact.")
GES (score-based, Gaussian) also returns a CPDAG: directed: ['X0->X2', 'X1->X2', 'X2->X3'] undirected: ['X0-X4']
Causal-sufficiency failure: A and B share a HIDDEN cause H (no true A-B edge). PC (not seeing H) infers a direct edge between A and B: True -> a phantom edge from unmeasured confounding. This is why discovery needs causal sufficiency (or FCI), and why a learned edge is a hypothesis to test, not a fact.
4. Real data — the Sachs protein-signaling network¶
Sachs et al. (Science, 2005) measured 11 phosphorylated proteins and phospholipids in thousands of individual human immune cells by flow cytometry, and used the data to reconstruct the cells' signaling network — the canonical real-world causal-discovery benchmark, because decades of molecular biology give it an established ground-truth DAG (17 directed edges: the Raf→Mek→Erk cascade, PKA/PKC as broad regulators, the Plcg→PIP3→PIP2 messenger chain).
We run PC on the observational subset (853 cells, log-transformed abundances) and score the recovered skeleton against the consensus network.
The result is better than "high precision, modest recall" suggests, and worse. Precision is perfect — not a single false edge — but that is a consequence of how little PC is willing to assert: 7 adjacencies out of 55 possible pairs. So we ask whether the missing edges are near misses, by sweeping the test level. They are not. Loosening the test twentyfold adds false edges and no true ones, which means the missed cascade is not hiding just past a significance threshold; it is absent from the observational distribution altogether. Discovery narrows the hypothesis space; it does not hand you the network, and tuning will not change that.
sachs=pd.read_csv("sachs.data.txt", sep="\t"); prot=list(sachs.columns)
Xs=np.log(sachs.values) # log abundances (standard for flow-cytometry)
print(f"Sachs single-cell data: {sachs.shape[0]} cells x {sachs.shape[1]} proteins {prot}")
# consensus ground-truth network (Sachs et al. 2005), 17 directed edges
truth=[("PKC","Raf"),("PKC","Mek"),("PKC","P38"),("PKC","Jnk"),("PKC","PKA"),
("PKA","Raf"),("PKA","Mek"),("PKA","Erk"),("PKA","Akt"),("PKA","P38"),("PKA","Jnk"),
("Raf","Mek"),("Mek","Erk"),("Erk","Akt"),
("Plcg","PIP2"),("Plcg","PIP3"),("PIP3","PIP2")]
truth_adj=set(frozenset(e) for e in truth)
cg=pc(Xs, 0.05, show_progress=False); G=cg.G.graph
found=set()
for i in range(len(prot)):
for j in range(i+1,len(prot)):
if G[i,j]!=0 or G[j,i]!=0: found.add(frozenset((prot[i],prot[j])))
tp=len(truth_adj & found); fp=len(found-truth_adj); fn=len(truth_adj-found)
print(f"\nPC vs the consensus network (skeleton / adjacency):")
print(f" true edges recovered (TP) = {tp}/{len(truth_adj)} false edges (FP) = {fp} missed (FN) = {fn}")
print(f" precision = {tp/(tp+fp):.2f} recall = {tp/(tp+fn):.2f}")
# network visualization
posS={"PKC":(0,3),"PKA":(2,3),"Raf":(1,2.1),"Mek":(1,1.1),"Erk":(1.3,0.15),"Akt":(2.3,0.15),
"Plcg":(3.7,3),"PIP3":(3.7,2),"PIP2":(3.7,1),"P38":(-0.4,1.4),"Jnk":(-0.4,0.35)}
def dnet(ax, edges_col, title):
for p,(x,y) in posS.items():
ax.scatter([x],[y],s=900,facecolor="white",edgecolor="k",zorder=3); ax.text(x,y,p,ha="center",va="center",fontsize=8,zorder=4)
for u,v,c,st in edges_col:
x1,y1=posS[u]; x2,y2=posS[v]; dx,dy=x2-x1,y2-y1; L=np.hypot(dx,dy); ux,uy=dx/L,dy/L; r=0.16
ax.plot([x1+ux*r,x2-ux*r],[y1+uy*r,y2-uy*r],color=c,lw=2,ls=st,zorder=2)
ax.set_title(title,fontsize=11); ax.axis("off"); ax.set_xlim(-0.9,4.2); ax.set_ylim(-0.3,3.4)
fig,ax=plt.subplots(1,2,figsize=(13,5))
dnet(ax[0],[(list(e)[0],list(e)[1],"k","-") for e in truth_adj],"Consensus network (17 edges, ground truth)")
col=[]
for e in found:
u,v=list(e); col.append((u,v,GREEN if e in truth_adj else RED,"-"))
for e in truth_adj-found:
u,v=list(e); col.append((u,v,GREY,":"))
dnet(ax[1],col,f"PC recovery: {tp} correct (green), {fp} false (red), {fn} missed (grey)")
plt.tight_layout(); plt.show()
npairs=len(prot)*(len(prot)-1)//2
print(f"Precision is not high here, it is perfect: {fp} false edges. That is worth explaining rather than")
print(f"praising, because the explanation is unflattering. PC asserted {len(found)} adjacencies out of {npairs} possible")
print(f"pairs, on data whose largest off-diagonal correlation is {np.abs(np.corrcoef(Xs.T)-np.eye(len(prot))).max():.2f}. It is not being")
print("accurate so much as being reluctant, and what it does commit to is safe.")
print()
print("So the question is whether the missing edges are near misses. Vary the test level:")
print()
print(f" {'alpha':>7} {'edges':>6} {'TP':>4} {'FP':>4} {'FN':>4} {'precision':>10} {'recall':>8}")
sweep=[]
for al in (0.001, 0.01, 0.05, 0.10, 0.20):
g=pc(Xs, al, show_progress=False).G.graph
f=set()
for i in range(len(prot)):
for j in range(i+1,len(prot)):
if g[i,j]!=0 or g[j,i]!=0: f.add(frozenset((prot[i],prot[j])))
t_=len(truth_adj & f); p_=len(f-truth_adj); m_=len(truth_adj-f)
sweep.append((al,len(f),t_,p_,m_))
print(f" {al:>7} {len(f):>6} {t_:>4} {p_:>4} {m_:>4} {t_/(t_+p_) if t_+p_ else float('nan'):>10.2f} {t_/(t_+m_):>8.2f}")
print()
tp_lo, tp_hi = sweep[2][2], sweep[-1][2]
print(f"The recall column does not move. Loosening the test twentyfold, from 0.01 to 0.20, adds")
print(f"{sweep[-1][1]-sweep[1][1]} edges and {sweep[-1][3]-sweep[1][3]} of them are false -- true positives go from {sweep[1][2]} to {tp_hi}. The edges PC misses")
print("are not sitting just past a significance boundary waiting for a more permissive threshold.")
print("They are invisible to conditional-independence testing on this data, and no choice of alpha")
print("recovers them. That is a statement about what the observational distribution contains, not")
print("about how the test was tuned.")
print()
Xraw=sachs.values.astype(float)
graw=pc(Xraw, 0.05, show_progress=False).G.graph
fraw=set()
for i in range(len(prot)):
for j in range(i+1,len(prot)):
if graw[i,j]!=0 or graw[j,i]!=0: fraw.add(frozenset((prot[i],prot[j])))
traw=len(truth_adj & fraw); praw=len(fraw-truth_adj)
print(f"One preprocessing note, since it was taken as standard rather than checked. On raw abundances")
print(f"instead of logs, PC finds {len(fraw)} edges: {traw} true, {praw} false, recall {traw/len(truth_adj):.2f} against {tp/len(truth_adj):.2f}. The log")
print("transform is conventional for flow cytometry and it costs a true edge here. Not a large effect,")
print("and worth knowing that the conventional choice was a choice.")
print()
print("Sachs et al. needed INTERVENTIONAL experiments -- perturbing each protein directly -- to recover")
print("the full directed network. That is the empirical version of this arc's through-line: the")
print("field's own benchmark for learning causal structure from observation was settled by")
print("intervening instead.")
Sachs single-cell data: 853 cells x 11 proteins ['Raf', 'Mek', 'Plcg', 'PIP2', 'PIP3', 'Erk', 'Akt', 'PKA', 'PKC', 'P38', 'Jnk'] PC vs the consensus network (skeleton / adjacency): true edges recovered (TP) = 7/17 false edges (FP) = 0 missed (FN) = 10 precision = 1.00 recall = 0.41
Precision is not high here, it is perfect: 0 false edges. That is worth explaining rather than
praising, because the explanation is unflattering. PC asserted 7 adjacencies out of 55 possible
pairs, on data whose largest off-diagonal correlation is 0.82. It is not being
accurate so much as being reluctant, and what it does commit to is safe.
So the question is whether the missing edges are near misses. Vary the test level:
alpha edges TP FP FN precision recall
0.001 6 6 0 11 1.00 0.35
0.01 6 6 0 11 1.00 0.35
0.05 7 7 0 10 1.00 0.41
0.1 7 7 0 10 1.00 0.41
0.2 10 7 3 10 0.70 0.41
The recall column does not move. Loosening the test twentyfold, from 0.01 to 0.20, adds
4 edges and 3 of them are false -- true positives go from 6 to 7. The edges PC misses
are not sitting just past a significance boundary waiting for a more permissive threshold.
They are invisible to conditional-independence testing on this data, and no choice of alpha
recovers them. That is a statement about what the observational distribution contains, not
about how the test was tuned.
One preprocessing note, since it was taken as standard rather than checked. On raw abundances
instead of logs, PC finds 8 edges: 8 true, 0 false, recall 0.47 against 0.41. The log transform is conventional for flow cytometry and it costs a true edge here. Not a large effect, and worth knowing that the conventional choice was a choice. Sachs et al. needed INTERVENTIONAL experiments -- perturbing each protein directly -- to recover the full directed network. That is the empirical version of this arc's through-line: the field's own benchmark for learning causal structure from observation was settled by intervening instead.
5. Summary¶
Causal discovery asks whether the DAG itself can be learned from data, and both the simulation and the real Sachs network set its promise and its limits:
- Constraint-based (PC) and score-based (GES) methods recover the skeleton and orient colliders (plus logical propagation), but return a CPDAG — the Markov equivalence class — leaving some edges undirected because Gaussian observational data cannot distinguish their direction.
- LiNGAM exploits non-Gaussianity to break those ties and recover the full DAG; on the simulation it correctly oriented the edge PC left ambiguous.
- On the real Sachs single-cell data, PC recovered the consensus network with perfect precision and 41% recall — every edge it proposed was genuine, because it proposed only 7 of 55 possible pairs. Sweeping the test level twentyfold added false edges and no true ones, so the missing cascade is not a tuning problem: it is not in the observational distribution. Sachs et al. needed interventions to reconstruct the full directed network.
- All of it rests on causal sufficiency (no hidden confounders — we showed a hidden common cause producing a phantom edge), faithfulness, and reliable independence testing.
The practical stance: treat causal discovery as a hypothesis generator, not a substitute for the identification designs and domain knowledge that make a causal claim credible — and note that the field's own benchmark was solved with interventions, echoing this arc's through-line that identification comes from design. Cross-links: this is the inverse of the DAGs & SCM notebook (subsection 8) — there the graph was given and we read off identification; here we learn only its equivalence class; Markov equivalence is the structural face of d-separation, and causal sufficiency is the discovery version of unconfoundedness; the interventions that complete the Sachs network are the molecular analogue of the randomized experiments in subsection 1. This completes the depth of the DAGs subsection. The R companion runs the same discovery with bnlearn (PC and hill-climbing), scoring against the Sachs ground truth.