Causal Inference VIII(b) — Causal Discovery (R companion)¶
bnlearn — constraint-based (PC) and score-based (hill-climbing) structure learning, on simulation and the Sachs network¶
bnlearn is R's standard structure-learning package. On the same Gaussian simulation as the Python notebook (true DAG $X_4\to X_0\to X_2$, $X_1\to X_2$, $X_2\to X_3$) it reproduces the central lesson — constraint-based pc.stable returns the CPDAG (leaving $X_0-X_4$ undirected), score-based hc returns a single DAG that may orient it wrongly — and it then applies both to the real Sachs (2005) protein-signaling data, scoring the recovered network against the established 17-edge ground truth. (LiNGAM lives in Python's causal-learn.)
1. Constraint-based PC — the CPDAG¶
pc.stable runs conditional-independence tests to recover the skeleton and orient v-structures. Its output is a CPDAG: the collider $X_0\to X_2\leftarrow X_1$ and $X_2\to X_3$ are directed, but $X_0-X_4$ is undirected — bnlearn represents an undirected edge as both X0->X4 and X4->X0.
suppressMessages(library(bnlearn))
d<-read.csv("discovery_gauss.csv")
pc<-pc.stable(d); ar<-arcs(pc)
cat("TRUE DAG: X4->X0, X0->X2, X1->X2, X2->X3\n\nPC (pc.stable) arcs:\n"); print(ar)
cat("\n-> X0-X4 appears in BOTH directions = undirected (unidentified); the rest are directed.\n")
pos<-list(X0=c(1,1),X1=c(2,0),X2=c(2,1),X3=c(3,1),X4=c(0,1))
drawdag<-function(ar,title){ plot(NA,xlim=c(-0.5,3.5),ylim=c(-0.4,1.5),axes=FALSE,xlab="",ylab="",main=title)
key<-apply(ar,1,function(r) paste(sort(r),collapse="_")); dup<-key %in% key[duplicated(key)]
for(k in seq_len(nrow(ar))){ a<-pos[[ar[k,1]]]; b<-pos[[ar[k,2]]]; u<-(b-a)/sqrt(sum((b-a)^2))*0.17
if(dup[k]){ if(ar[k,1]<ar[k,2]) lines(c(a[1]+u[1],b[1]-u[1]),c(a[2]+u[2],b[2]-u[2]),col="#c53030",lwd=3,lty=2) }
else arrows(a[1]+u[1],a[2]+u[2],b[1]-u[1],b[2]-u[2],length=0.12,lwd=2) }
for(nm in names(pos)){ p<-pos[[nm]]; points(p[1],p[2],pch=21,bg="white",cex=5); text(p[1],p[2],nm) } }
options(repr.plot.width=7, repr.plot.height=4)
drawdag(ar, "PC learns the CPDAG (X0-X4 undirected, dashed red)")
Warning message: "package 'bnlearn' was built under R version 4.6.1"
TRUE DAG: X4->X0, X0->X2, X1->X2, X2->X3 PC (pc.stable) arcs:
from to [1,] "X0" "X2" [2,] "X0" "X4" [3,] "X1" "X2" [4,] "X2" "X3" [5,] "X4" "X0"
-> X0-X4 appears in BOTH directions = undirected (unidentified); the rest are directed.
2. Score-based hill-climbing — a single DAG from the class¶
hc greedily maximizes BIC over DAG space and returns one DAG. Forced to output a fully-directed graph, it picks a direction for the ambiguous $X_0-X_4$ edge — and, since both directions are score-equivalent, it may pick the wrong one. The CPDAG is the more honest object.
hc_<-hc(d); arh<-arcs(hc_)
cat("hill-climbing (hc, score-based BIC) arcs (a single DAG):\n"); print(arh)
x04<-arh[apply(arh,1,function(r) all(sort(r)==c("X0","X4"))),,drop=FALSE]
cat(sprintf("\nhc oriented the X0-X4 edge as %s->%s", x04[1,1], x04[1,2]))
cat(if(x04[1,1]=="X0") " -- the WRONG direction (true is X4->X0); score-equivalent, so hc guessed.\n" else " (correct here, but not guaranteed).\n")
options(repr.plot.width=7, repr.plot.height=4)
drawdag(arh, "hc returns one DAG from the equivalence class")
hill-climbing (hc, score-based BIC) arcs (a single DAG):
from to [1,] "X2" "X3" [2,] "X0" "X4" [3,] "X0" "X2" [4,] "X1" "X2"
hc oriented the X0-X4 edge as X0->X4
-- the WRONG direction (true is X4->X0); score-equivalent, so hc guessed.
3. Real data — the Sachs protein-signaling network¶
Sachs et al. (2005) measured 11 signaling proteins in thousands of single immune cells; molecular biology gives the network an established 17-edge ground-truth DAG. We run pc.stable (partial-correlation tests) and hc (Gaussian BIC) on the observational log-abundances and score the recovered skeleton against that consensus. As in Python, the result is high precision, modest recall — observational data alone recovers a genuine subnetwork but misses much of the cascade that Sachs resolved with interventions.
s<-read.csv("sachs.data.txt", sep="\t"); s<-as.data.frame(log(s)); prot<-names(s)
cat(sprintf("Sachs single-cell data: %d cells x %d proteins\n", nrow(s), ncol(s)))
truth<-list(c("PKC","Raf"),c("PKC","Mek"),c("PKC","P38"),c("PKC","Jnk"),c("PKC","PKA"),
c("PKA","Raf"),c("PKA","Mek"),c("PKA","Erk"),c("PKA","Akt"),c("PKA","P38"),c("PKA","Jnk"),
c("Raf","Mek"),c("Mek","Erk"),c("Erk","Akt"),
c("Plcg","PIP2"),c("Plcg","PIP3"),c("PIP3","PIP2"))
tkey<-sapply(truth,function(e) paste(sort(e),collapse="_"))
score<-function(ar,name){
sk<-unique(apply(ar,1,function(r) paste(sort(r),collapse="_")))
tp<-sum(sk %in% tkey); fp<-sum(!(sk %in% tkey)); fn<-sum(!(tkey %in% sk))
cat(sprintf(" %-10s TP=%2d/%d FP=%2d FN=%2d precision=%.2f recall=%.2f\n", name, tp, length(tkey), fp, fn, tp/(tp+fp), tp/(tp+fn)))
invisible(sk) }
pcS<-pc.stable(s, test="cor"); hcS<-hc(s, score="bic-g")
cat("Recovery vs the 17-edge consensus network (skeleton):\n")
sk_pc<-score(arcs(pcS),"pc.stable"); sk_hc<-score(arcs(hcS),"hc")
# network plot of pc.stable recovery vs truth
posS<-list(PKC=c(0,3),PKA=c(2,3),Raf=c(1,2.1),Mek=c(1,1.1),Erk=c(1.3,0.15),Akt=c(2.3,0.15),
Plcg=c(3.7,3),PIP3=c(3.7,2),PIP2=c(3.7,1),P38=c(-0.4,1.4),Jnk=c(-0.4,0.35))
dnet<-function(sk,title){ plot(NA,xlim=c(-0.9,4.2),ylim=c(-0.3,3.4),axes=FALSE,xlab="",ylab="",main=title)
alledges<-union(sk,tkey)
for(k in alledges){ pr<-strsplit(k,"_")[[1]]; a<-posS[[pr[1]]]; b<-posS[[pr[2]]]; u<-(b-a)/sqrt(sum((b-a)^2))*0.16
col<-if(k %in% sk && k %in% tkey) "#2f855a" else if(k %in% sk) "#c53030" else "#a0aec0"
lty<-if(k %in% sk) 1 else 3
lines(c(a[1]+u[1],b[1]-u[1]),c(a[2]+u[2],b[2]-u[2]),col=col,lwd=2,lty=lty) }
for(nm in names(posS)){ p<-posS[[nm]]; points(p[1],p[2],pch=21,bg="white",cex=3.6); text(p[1],p[2],nm,cex=.7) } }
options(repr.plot.width=8, repr.plot.height=5.5)
dnet(sk_pc, "pc.stable on Sachs: green=correct, red=false, grey dashed=missed")
cat("\nPrecision is perfect, not merely high: both algorithms proposed 7 adjacencies out of 55 possible pairs and every one is real.
That is reluctance rather than accuracy -- see the Python notebook, where sweeping the test level twentyfold adds false
edges and no true ones. Observational discovery narrows the hypothesis space; Sachs needed interventions for the rest.\n")
Sachs single-cell data: 853 cells x 11 proteins
Recovery vs the 17-edge consensus network (skeleton):
pc.stable TP= 7/17 FP= 0 FN=10 precision=1.00 recall=0.41
hc TP= 7/17 FP= 0 FN=10 precision=1.00 recall=0.41
Precision is perfect, not merely high: both algorithms proposed 7 adjacencies out of 55 possible pairs and every one is real. That is reluctance rather than accuracy -- see the Python notebook, where sweeping the test level twentyfold adds false edges and no true ones. Observational discovery narrows the hypothesis space; Sachs needed interventions for the rest.
4. Summary¶
bnlearn reproduced the structure-learning results across simulation and real data: pc.stable returned the CPDAG with $X_0-X_4$ undirected, hc returned a single DAG that had to guess that edge (and can get it wrong), and on the real Sachs network both recovered a high-precision, modest-recall subnetwork of the 17-edge consensus. The direction PC leaves open is the one LiNGAM (Python causal-learn) orients via non-Gaussianity; the edges observational discovery misses are the ones Sachs resolved with interventions.
The lesson holds across languages and methods: observational discovery yields an equivalence class and a partial network, not a unique causal graph, and every learned structure rests on causal sufficiency and faithfulness. Use it to generate hypotheses, then confirm with domain knowledge, experiments, or the identification designs of this arc. Cross-links: the inverse of the DAGs & SCM notebook; Markov equivalence is the structural counterpart of d-separation; causal sufficiency is the discovery form of unconfoundedness; the Sachs interventions echo the randomized experiments of subsection 1. This completes the depth of the DAGs subsection.