Scalable Bayesian structure learning of directed acyclic graphs via Laplace approximation, with an application to breast cancer gene expression networks
Structure learning of directed acyclic graphs (DAGs) from observational data is a foundational task in causal discovery and is widely used to infer regulatory networks from medical and genomic measurements. The Bayesian formulation quantifies model uncertainty and admits prior biological knowledge, but its practical use has been hampered by the super-exponential growth of the DAG space and by the intractability of the node-marginal likelihood under flexible, non-conjugate priors. Existing closed-form solutions are largely confined to the conjugate Normal--Inverse-Gamma prior. We develop a Laplace-approximated Bayesian scoring function for the non-conjugate Normal--Gamma prior on the modified Cholesky parameterisation of the precision matrix, embed it in a Metropolis--Hastings sampler over DAGs, and couple the latent Gaussian network to a binary clinical outcome through a probit link. We show that the node-marginal integral is of generalised inverse-Gaussian form, so that its exact value is a modified Bessel function of the second kind and the proposed scoring function is its leading large-argument asymptotic; the posterior of each conditional variance is likewise generalised inverse-Gaussian and is sampled exactly. In simulation, the proposed prior improves on the conjugate baseline and on the PC, greedy-equivalence-search, NOTEARS, and DAGMA benchmarks at sample sizes typical of clinical cohorts. On two real datasets, the Sachs protein-signalling network, scored against its validated consensus graph, and the Wisconsin Diagnostic Breast Cancer data, the method recovers known structure and, through the DAG-probit extension, predicts malignancy from nuclear morphometry with a cross-validated ROC-AUC of $0.94$ using a sparse, interpretable set of direct predictors.