Bayesian Causal Discovery - BGe Score, Interventions and Sachs Data

Summary

Concept note (no raw source). For continuous data, Bayesian DAG learning uses the BGe score (Normal-Wishart marginal likelihood; score-equivalent, so Markov-equivalent DAGs tie). Observational data identify only a Markov equivalence class; interventional data (Sachs et al. 2005: 7,466 cells over 9 conditions) break ties by removing the intervened node’s local factor. Posterior averaging over DAGs uses order-based sums (see Friedman 2000 - Closed Form for a Known Ordering) or dynamic programming over orders. Used in PyMC example E20.

Overview

Example list: https://redam94.github.io/pymc-challenges/examples.html, E20 (linear-Gaussian SEMs, Markov equivalence through LOO and pm.do, NOTEARS acyclicity as pm.Potential and why it is multimodal and scale-dependent, fixed-order horseshoe DAGs, exact BGe + dynamic-programming averaging over all orders). E19 is the undirected sibling (Bayesian Gaussian Graphical Models - G-Wishart and Neighbourhood Selection).

Main Content

Linear-Gaussian structural equation model on a DAG

$x_i={\sum }_{j\in \text{Pa}(i)}{w}_{ji}{x}_j+{ϵ}_i$, ${ϵ}_i\sim N(0,{\sigma }_i^2)$ independent; $W$ is the weighted adjacency matrix with acyclic support. Observational distribution $N(0,(I-W{)}^{-\top }{\Sigma }_{ϵ}(I-W{)}^{-1})$; different DAGs in one Markov equivalence class give the same distribution (e.g. $X\to Y$ and $Y\to X$).

BGe score (Geiger & Heckerman 1994, 2002)

With a Normal-Wishart prior ($\mu \mid \Omega \sim N({\mu }_0,({\alpha }_{\mu }\Omega {)}^{-1})$, $\Omega \sim W({\alpha }_w,T)$) over $m$ variables, the marginal likelihood of $n$ complete observations of a subset of $k$ variables is

$$p(D_k)=(2\pi {)}^{-nk/2}(\frac{{\alpha }_{\mu }}{{\alpha }_{\mu }+n}{)}^{k/2}\frac{c(k,{\alpha }_w^{(k)})}{c(k,{\alpha }_w^{(k)}+n)}|T_k{|}^{{\alpha }_w^{(k)}/2}|R_k{|}^{-({\alpha }_w^{(k)}+n)/2},$$

$R_k=T_k+S_n+\frac{{\alpha }_{\mu }n}{{\alpha }_{\mu }+n}({\mu }_0-\bar{x})({\mu }_0-\bar{x}{)}^{\top }$, $c(k,\alpha )=[2^{\alpha k/2}{\pi }^{k(k-1)/4}{\prod }_{i=1}^k\Gamma (\frac{\alpha +1-i}{2}){]}^{-1}$, ${\alpha }_w^{(k)}={\alpha }_w-m+k$. Local score $\text{score}(X_i,U\mid D)=p(D_{\{i\}\cup U})/p(D_U)$; $P(D\mid G)={\prod }_i\text{score}(X_i,{\text{Pa}}_G(X_i)\mid D)$, as in Friedman 2000 - Bayesian Structure Learning Framework. BGe is score equivalent: Markov-equivalent DAGs have equal scores.

Interventional data (Cooper & Yoo 1999; Sachs et al. 2005)

For samples where node $i$ is experimentally set (perfect intervention), $x_i$ is not generated by its parents, so that sample is excluded from node $i$‘s local score (its parent mechanism is not observed), while its value still acts as a parent for children. In PyMC, pm.do / per-sample masks implement the same idea: the likelihood term for $x_i$ in intervened rows is dropped. Interventions make edge directions identifiable beyond the Markov class.

Sachs et al. (2005): flow cytometry on 11 phosphoproteins/phospholipids in primary T cells under multiple perturbations; the commonly used benchmark has a 20-edge consensus network and (per the example list) 7,466 cells in 9 intervention conditions.

Computation: (a) order-MCMC / exact sums over DAGs given an order (Friedman & Koller); (b) exact averaging over all orders by dynamic programming over subsets (Koivisto & Sood 2004; cost about $O(n{2}^n)$ local-score lookups with an in-degree bound), practical to roughly $n\lesssim 15$-20; (c) continuous relaxations: NOTEARS $h(W)=\text{tr}{e}^{W\circ W}-d=0$ (NOTEARS - Overview), which in a Bayesian model becomes pm.Potential(-rho*h(W)): soft, scale-dependent and multimodal.

Examples

Score equivalence check ( $X\to Y$ vs $Y\to X$)

$n=200$, $x\sim N(0,1)$, $y=0.8x+ϵ$. Both orientations get BGe log-score $-580.5959$ (identical), the empty graph $-615.7926$: data support an edge but not its direction.

import numpy as np
from scipy.special import gammaln
def logc(k, a): return -(a*k/2*np.log(2) + k*(k-1)/4*np.log(np.pi) + sum(gammaln((a+1-i)/2) for i in range(1, k+1)))
def bge_node(D, i, pa, am=1.0):
    n, m = D.shape; aw_full = m + am + 1
    def lm(cols):
        X = D[:, cols]; k = len(cols); aw = aw_full - m + k
        T = (am*(aw-k-1)/(am+1))*np.eye(k); xb = X.mean(0); Sn = (X-xb).T @ (X-xb)
        R = T + Sn + (am*n/(am+n))*np.outer(xb, xb)
        return (-n*k/2*np.log(2*np.pi) + k/2*np.log(am/(am+n)) + logc(k, aw) - logc(k, aw+n)
                + aw/2*np.linalg.slogdet(T)[1] - (aw+n)/2*np.linalg.slogdet(R)[1])
    return lm([i] + list(pa)) - (lm(list(pa)) if pa else 0.0)
rng = np.random.default_rng(0); x = rng.normal(size=200); y = 0.8*x + rng.normal(size=200); D = np.c_[x, y]
print(bge_node(D,0,[]) + bge_node(D,1,[0]), bge_node(D,1,[]) + bge_node(D,0,[1]))

(Prior: ${\mu }_0=0$, $T=\frac{{\alpha }_{\mu }({\alpha }_w-m-1)}{{\alpha }_{\mu }+1}I$, ${\alpha }_{\mu }=1$, ${\alpha }_w=m+{\alpha }_{\mu }+1$, standard defaults.)

Marketing relevance: causal structure among marketing variables (spend, brand search, awareness, sales) is only partially identifiable from observational time series; experiments (geo holdouts, lift tests) play the role of interventions, and a posterior over edges communicates that honestly. Treat as an analogy; time-series dependence is not handled by the iid BGe model.

Connections

Causal Discovery - Overview, PC Algorithm and Constraint-Based Discovery, GES and Score-Based Discovery, Markov and Faithfulness Assumptions, Functional Causal Models (LiNGAM, ANM), Joint Priors and Covariance Matrices.

See Also

Friedman 2000 - Being Bayesian about Network Structure