Home | About | Journals | Submit | Contact Us | Français |

**|**Europe PMC Author Manuscripts**|**PMC5358773

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Notation and Model
- 3. A Gibbs Sampler for Structure Learning
- 4. Computational Aspects
- 5. Results
- 6. Discussion
- Supplementary Material
- References

Authors

Related links

J Mach Learn Res. Author manuscript; available in PMC 2017 March 20.

Published in final edited form as:

J Mach Learn Res. 2016 April; 17(30): 1–39.

PMCID: PMC5358773

EMSID: EMS67582

Peter Spirtes, Editor

Robert J. B. Goudie, Medical Research Council Biostatistics Unit Cambridge CB2 0SR, UK;

We propose a Gibbs sampler for structure learning in directed acyclic graph (DAG) models. The standard Markov chain Monte Carlo algorithms used for learning DAGs are random-walk Metropolis-Hastings samplers. These samplers are guaranteed to converge asymptotically but often mix slowly when exploring the large graph spaces that arise in structure learning. In each step, the sampler we propose draws entire sets of parents for multiple nodes from the appropriate conditional distribution. This provides an efficient way to make large moves in graph space, permitting faster mixing whilst retaining asymptotic guarantees of convergence. The conditional distribution is related to variable selection with candidate parents playing the role of covariates or inputs. We empirically examine the performance of the sampler using several simulated and real data examples. The proposed method gives robust results in diverse settings, outperforming several existing Bayesian and frequentist methods. In addition, our empirical results shed some light on the relative merits of Bayesian and constraint-based methods for structure learning.

We consider structure learning for graphical models based on directed acyclic graphs (DAGs). The basic structure learning task can be stated simply: given data **X** on *p* variables, assumed to be generated from a graphical model based on an unknown DAG *G*, the goal is to make inferences concerning the edge set of *G* (the vertex set is treated as fixed and known). Structure learning appears in a variety of applications (for an overview see Korb and Nicholson, 2011) and is a key subtask in many analyses involving graphical models, including causal inference.

In a Bayesian framework, structure learning is based on the posterior distribution *P*(*G* | **X**) (Koller and Friedman, 2009). The domain of the distribution is the space $\mathcal{G}$ of all DAGs with *p* vertices. The size of the space grows super-exponentially with *p*, precluding exhaustive enumeration for all but the smallest problems. Markov chain Monte Carlo (MCMC) methods are widely used to sample DAGs and thereby approximate the posterior *P*(*G* | **X**). Available methods include MC^{3} (Madigan and York, 1995) and related samplers (Giudici and Castelo, 2003; Grzegorczyk and Husmeier, 2008) and algorithms that sample from the space of total orders (Friedman and Koller, 2003). Order-space sampling entails some restrictions on graph priors (Eaton and Murphy, 2007; Ellis and Wong, 2008), because the number of total orders with which a DAG is consistent is not constant. Order-space approaches have more recently led to exact methods based on dynamic programming (Koivisto and Sood, 2004; Parviainen and Koivisto, 2009; Tian and He, 2009; Tamada et al., 2011; Parviainen and Koivisto, 2013; Nikolova et al., 2013). These are currently feasible only in small domain settings (typically *p* < 20 but *p* < 30 is feasible on large cluster computers) and usually share the same restrictions on graph priors as order-space samplers. Frequentist constraint-based methods such as the PC-algorithm (Spirtes et al., 2000; Colombo and Maathuis, 2014) and the approach of Xie and Geng (2008) are an alternative. These methods make firm decisions about the DAG structure via a series of tests of conditional independence.

In this paper we propose a Gibbs sampler for structure learning of DAGs that ameliorates key deficiencies in existing samplers. Random-walk Metropolis-Hastings algorithms currently used for structure learning propose ‘small’ changes at each iteration, which can leave the algorithms unable to escape local modes. The Gibbs sampler proposed here considers the parents of a set of nodes as a single component (a so-called ‘block’), sampling an entire parent set for each node in the block in one step. These ‘large’ moves are sampled from the appropriate conditional posterior distribution. This enables the sampler to efficiently locate and explore areas of significant posterior mass. The method is based on the simple heuristic that the parents of a node are similar to the covariates or inputs in variable selection, with the node as the output variable, but accounts for acyclicity exactly so that the equilibrium distribution is indeed the correct posterior distribution over DAGs. The sampler does not impose restrictions on priors or graph space beyond those (maximum in-degree, modular prior) common to most samplers for structure learning (e.g. Friedman and Koller, 2003; Ellis and Wong, 2008; Grzegorczyk and Husmeier, 2008). The maximum in-degree restriction formally precludes large-sample consistency in the general case, but facilitates effective and robust inference by reducing the size of the model space (see Discussion).

Let *G* denote a DAG with vertex set *V* (*G*) = {1, … , *p*}, and directed edge set *E*(*G*) *V* (*G*) × *V* (*G*); often we will refer to vertex and edge sets simply as *V* and *E* respectively, leaving *G* implicit. The binary adjacency matrix corresponding to *G* is denoted *A ^{G}*, with entries specified as ${A}_{uv}^{G}=1\iff (u,v)\in E(G)$ and diagonal entries equal to zero. For inference,

The proposed Gibbs sampler changes entire parent sets *en masse*, not just for one node but for a set of nodes, and so it will often be convenient and natural to specify a DAG *G* in terms of the parents of each node. Let ${\text{pa}}_{v}^{G}=\{u:(u,v)\in E(G)\}$ denote the set of nodes that are parents of node *v* in *G*. A tuple of *p* parent sets $\langle {\text{pa}}_{1}={\text{pa}}_{1}^{G},\dots ,{\text{pa}}_{p}={\text{pa}}_{p}^{G}\rangle $ fully specifies the edge set *E*, since $u\in {\text{pa}}_{v}^{G}\iff \left(u,v\right)\in E(G).$ Thus, structure learning can be viewed as inference of parent sets pa_{v}. The parent set pa_{v} takes values in the power set $\mathcal{P}$(*V* \ {*v*}) of nodes excluding node *v*, subject to acyclicity. We will usually suppress the labels, and simply write $\langle {\text{pa}}_{1}^{G},\dots ,{\text{pa}}_{p}^{G}\rangle $ when specifying parent sets. Since ${\text{pa}}_{v}^{G}$ is a subset of the variable indices, we can use ${\mathbf{\text{X}}}_{{\text{pa}}_{v}^{G}}$ to denote the corresponding data submatrix.

We denote the tuple of parent sets of a set of nodes *Z* = {*z*_{1}, … , *z _{s}*}

Illustration of notation. A graph *G*, with vertex set *V* = {1, 2, 3, 4}, edge set *E* = {(1, 2), (1, 3), (2, 3), (4, 3)} and adjacency matrix *A*^{G} as shown, can be specified using parent sets as *G* = pa_{1} = , {1}, pa_{2} = {1}, pa_{3} ={1,2,4}, pa **...**

Our statistical formulation for Bayesian networks is standard. We briefly summarize the main points and refer the reader to the references below for details. Given a DAG *G*, the joint distribution of *X* is $p\left(X|G,\left\{{\theta}_{v}\right\}\right)={\displaystyle {\prod}_{v\in V}p\left({X}_{v}|{X}_{{\text{pa}}_{v}^{G}},{\theta}_{v}\right),}$ i.e. the joint distribution factors over nodes, and each node is conditioned on its parents in *G*, parameterized by *θ _{v}*. For structural inference, interest focuses on the posterior distribution

$$\begin{array}{l}P\left(G|\mathbf{\text{X}}\right)=P({\text{pa}}_{1}={\text{pa}}_{1}^{G},\dots ,{\text{pa}}_{p}={\text{pa}}_{p}^{G}|\mathbf{\text{X}})\\ \phantom{\rule{3.2em}{0ex}}\propto {\displaystyle \prod _{v\in V}p({\mathbf{\text{X}}}_{v}|{\mathbf{\text{X}}}_{\text{p}{\text{a}}_{v}^{G}}){\pi}_{v}}({\text{pa}}_{v}^{G}),\end{array}$$

where $p\left({\mathbf{\text{X}}}_{v}|{\mathbf{\text{X}}}_{{\text{pa}}_{v}^{G}}\right)$ is the marginal likelihood for node *v* given the graph $G=\langle {\text{pa}}_{1}^{G},\dots ,{\text{pa}}_{p}^{G}\rangle .$ The distribution *P* (*G* | **X**) is the target distribution for our sampler.

In this section we provide a high-level description of the Gibbs sampler for structure learning. We first recall the standard random-walk Metropolis-Hastings sampler (known as MC^{3}) and then describe a naïve Gibbs sampler, which offers no gains over MC^{3}, but prepares the ground for the introduction of ideas from the Gibbs sampling literature. Specifically we discuss a strategy known as blocking that can improve mixing in Gibbs samplers and show how to use blocking for structure learning of DAGs. Several points relating to computation are central to developing a practical Gibbs sampler in this setting, but for clarity of exposition we defer discussion of computational aspects to Section 4.

The standard sampler for structure learning of DAGs is MC^{3} (Madigan and York, 1995). This is a classical Metropolis-Hastings sampler with proposal *G*′ drawn uniformly at random from the set neigh(*G*) of DAGs that differ from the current DAG *G* by the addition or removal of a single edge. The acceptance probability is min(1, *r*(*G*′, *G*)), where

$$r({G}^{\prime},G)=\mathrm{min}\left\{1,\frac{P({G}^{\prime}|\mathbf{\text{X}})|\text{neigh}({G}^{\prime}){|}^{-1}}{P(G|\mathbf{\text{X}})|\text{neigh}(G){|}^{-1}}\right\}.$$

Variants of MC^{3} include single-edge direction reversal proposals (Giudici and Castelo, 2003).

Constructing a Gibbs sampler that is analogous to MC^{3} is straightforward. The posterior distribution on DAGs is a joint distribution over the off-diagonal entries in the adjacency matrix *A ^{G}*, i.e. over the

A simple Gibbs sampler works in a similar way. At each step, a move from graph *G* to a new graph *G*′ is chosen by sampling from the conditional distribution of ${A}_{uv}^{{G}^{\prime}}$ (for some *u* ≠ *v*) given the rest of the graph. If (*u, v*) *E*(*G*), define the graph *G*^{(u,v)} to be identical to *G*, and *G*^{−(u,v)} to be the graph that differs from *G* only in lacking an edge from *u* to *v*; conversely, if (*u, v*) *E*(*G*), define *G*^{−(u,v)} to be identical to *G*, and *G*^{(u,v)} to be the graph that differs from *G* only in including an edge from *u* to *v*. If *G*^{(u,v)} is cyclic, *G*^{−(u,v)} is sampled as the new graph *G*′ with probability 1. If *G*^{(u,v)} is acyclic, the conditional distribution of ${A}_{uv}^{{G}^{\prime}}$ is Bernoulli, i.e.

$$P({A}_{uv}^{{G}^{\prime}}=a\phantom{\rule{0.1em}{0ex}}|\phantom{\rule{0.1em}{0ex}}{A}_{-uv}^{G})=\{\begin{array}{cc}\frac{P({G}^{-(u,v)}|\mathbf{\text{X}})}{P({G}^{-(u,v)}|\mathbf{\text{X}})+P({G}^{(u,v)}|\mathbf{\text{X}})}& a=0,\\ \frac{P({G}^{(u,v)}|\mathbf{\text{X}})}{P({G}^{-(u,v)}|\mathbf{\text{X}})+P({G}^{(u,v)}|\mathbf{\text{X}})}& a=1.\end{array}$$

The choice of *u* and *v* can either be made sequentially (systematically) or randomly (Roberts and Sahu, 1997); in this paper, random-scan Gibbs samplers are used throughout. We prove convergence of the Gibbs sampler to the target distribution in Appendix A.

The mixing of Metropolis-Hastings samplers depends upon the proposal distribution, which for convenience is often chosen as a random-walk. In contrast, Gibbs samplers make moves according to conditional distributions that reflect local structure of the target distribution. Nonetheless, Gibbs sampling is not always efficient. In particular, correlation between the components being sampled can lead to inefficiency. To see this, consider a Gibbs sampler for a multivariate continuous distribution with highly correlated components. At each step, a single component is sampled according to its conditional distribution, but since it is strongly correlated with other component(s), the conditional is concentrated on only a small part of the support. This means the sampler is likely to make only small moves. Analogous issues arise with discrete distributions.

For graphical models based on DAGs, there may be strong dependence between the edge indicators ${A}_{uv}^{G}$, particularly for the collections of RVs corresponding to parent sets. For example, there may be RVs *X _{u}* and

We address the deficiencies of the naïve Gibbs sampler by grouping a number of the components together and sampling from their joint conditional distribution. In Gibbs sampling, this is known as *blocking*. Sampling from such a joint conditional can ameliorate difficulties caused by correlations between components because the joint conditional naturally accounts for the correlation structure. In the multivariate normal case, Roberts and Sahu (1997) have shown that blocking improves convergence for random-scan Gibbs sampling.

In principle, any group of components can be taken as a block, but for the algorithm to be useful in practice, sampling from a block’s joint conditional distribution must be feasible, and ideally simple. The blocks that we consider correspond to groups (specifically tuples) of parent sets, so that the parent sets of several nodes are considered simultaneously. It is natural that each block is a tuple of parent sets because, as described in Section 2, we can specify a graph *G* using a tuple $\langle {\text{pa}}_{1}^{G},\dots ,{\text{pa}}_{p}^{G}\rangle $ of parent sets. It is therefore convenient to describe the sampler using this alternative graph specification, in which $G=\langle {\text{pa}}_{1}^{G},\dots ,{\text{pa}}_{p}^{G}\rangle .$

We denote the set of *q* nodes whose parent sets will be sampled together as a block by *W* = {*w*_{1}, … , *w _{q}*}

To construct a Gibbs sampler using these blocks, we need to find the conditional posterior distribution on the tuple pa_{W} of parent sets, given that the remaining tuple pa_{−W} of parent sets is set to ${\text{pa}}_{-W}^{G}.$ The conditional distribution depends on whether the graph *G*′ formed using the proposed parent sets is acyclic. We therefore introduce ${P\text{a}}_{W}^{G}$ to denote the set of permissible tuples, i.e. tuples ${\text{pa}}_{W}^{{G}^{\prime}}$ of parent sets such that ${G}^{\prime}=\langle {\text{pa}}_{W}^{{G}^{\prime}},{\text{pa}}_{-W}^{G}\rangle $ is acyclic. For tuples ${\text{pa}}_{W}^{{G}^{\prime}}\notin P{\text{a}}_{W}^{G},$ the conditional probability is 0. For ${\text{pa}}_{W}^{{G}^{\prime}}\in P{\text{a}}_{W}^{G},$ the conditional probability is

$$P({\text{pa}}_{W}={\text{pa}}_{W}^{{G}^{\prime}}|{\text{pa}}_{-W}={\text{pa}}_{-W}^{G},\mathbf{\text{X}})=\frac{P(\langle {\text{pa}}_{W}^{{G}^{\prime}},{\text{pa}}_{-W}^{G}\rangle |\mathbf{\text{X}})}{P({\text{pa}}_{-W}^{G}|\mathbf{\text{X}})}\phantom{\rule{0ex}{0ex}}\phantom{\rule{14em}{0ex}}=\frac{P({G}^{\prime}|\mathbf{\text{X}})}{{\displaystyle {\sum}_{{\text{pa}}_{W}^{{G}^{\u2033}}\in P{\text{a}}_{W}^{G}}P({\text{pa}}_{W}^{{G}^{\u2033}},{\text{pa}}_{-W}^{G}|\mathbf{\text{X}})}}.$$

(1)

This is the conditional distribution needed to specify the blocked Gibbs sampler; Algorithm 1 outlines the procedure at a high level, with the set *W* of nodes chosen uniformly at random at each step. Asymptotic convergence follows from the argument in Appendix A (in fact the requirements on the graph prior will be weaker than in the naïve case).

Initialise starting point ${G}_{0}=\langle {\text{pa}}_{1}^{{G}_{0}},\dots ,{\text{pa}}_{p}^{{G}_{0}}\rangle $

**for**
*t* in 1 to *N*
**do**

Sample *q* nodes uniformly at random from *V*, and call this set of nodes *W*

Sample ${\text{pa}}_{W}^{{G}^{\prime}}\phantom{\rule{0.2em}{0ex}}\text{from}\hspace{0.17em}P({\text{pa}}_{W}={\text{pa}}_{W}^{{G}^{\prime}}|{\text{pa}}_{-W}^{{G}_{t-1}},\text{X})$

Set ${G}_{t}\leftarrow G=\langle {\text{pa}}_{W}^{{G}^{\prime}},{\text{pa}}_{-W}^{{G}_{t-1}}\rangle $

**end for**

To reduce the computational costs of structure learning it is common to set a maximum in-degree (e.g. Friedman and Koller, 2003; Grzegorczyk and Husmeier, 2008). We set a maximum in-degree *κ* = 3 in all empirical examples, except where stated otherwise (Section 5). This facilitates sampling from the conditional distribution in Equation 1 by reducing the computational cost of evaluating the normalising constant. We set the block size *q* = *|W|* = 3. While *q* can be chosen freely in principle, evaluating the conditional distribution in Equation 1 becomes unmanageable when *q* is large. Thus, both *κ* and *q* act as tuning parameters (and are in addition to any hyper-parameters in the Bayesian formulation).

It is interesting to note that when *q* = 1 (and thus *W* = *w* for some *w* *V*) if no choice of parent set induces a cycle then $P{\text{a}}_{W}^{G}=\mathcal{P}(V\backslash \{w\}),$ the power set of the remaining nodes. In this case, the conditional distribution in Equation 1 can be viewed as the posterior distribution of a variable selection problem with response or output variable *w*, and the other variables as covariates or inputs. If the addition of particular nodes introduces a cycle, we have a constrained problem. In particular, suppose adding any of the nodes *Z* *V* as parents of node *w* would create a cycle. Then $P{\text{a}}_{W}^{G}=\mathcal{P}(V\backslash \{w\}\cup Z))$ and the conditional distribution in Equation 1 is a constrained variable selection problem in which the nodes *Z* are excluded from the set of candidate inputs.

Designing a computationally efficient sampler is not straightforward in this setting. To see why, note that to sample from the conditional in Equation 1 we need to be able to identify the set $P{\text{a}}_{W}^{G}$ of tuples of parent sets that is permissible (in the sense of maintaining acyclicity). The cardinality of this set is typically large, and the interdependence between the parent sets of each node in *W* makes decomposing the problem into subproblems non-trivial. A simple but naïve approach would list all possible parent sets for each node in *W* and check each such combination for cyclicity, but this approach is slow and cumbersome.

We propose a partitioning scheme that leads to a two-stage sampler. The key idea is to choose the partition of $P{\text{a}}_{W}^{G}$ so that, conditional on a component of the partition, the parents of each node are independent. This enables an efficient two-stage sampling method that we describe in Section 4.3. Acyclicity constraints are met by an efficient dynamic algorithm.

We partition the set $P{\text{a}}_{W}^{G}$ of permissible tuples of parent sets as $P{\text{a}}_{W}^{G}=\left\{P{\text{a}}_{W}^{G,{H}_{1}},\dots ,P{\text{a}}_{W}^{G,{H}_{\eta}}\right\}.$ Each component $P{\text{a}}_{W}^{G,{H}_{h}}$ is a set of permissible tuples of parent sets for nodes in *W*. It is convenient to label the partition components using secondary (unrelated to *G*) DAGs *H _{h}* , where is the space of DAGs with

- The reduced graph $\overline{G}=\langle \text{p}{\text{a}}_{1}^{\overline{G}},\dots ,\text{p}{\text{a}}_{p}^{\overline{G}}\rangle ,$ which is a function of the graph
*G*, and is identical to*G*except that edges directed into nodes in*W*are removed.$${\text{pa}}_{w}^{\overline{G}}=\{\begin{array}{ll}\varnothing \hfill & w\in W,\hfill \\ {\text{pa}}_{w}^{G}\hfill & w\notin W\hfill \end{array}$$ - The (reflexive) transitive closure, which is the directed graph T
^{G}on nodes*V*with edges*E*^{T}, where (*u*,*v*)*E*^{T}if and only if*u*=*v*or a (directed) path from*u*to*v*exists in*G*(i.e. there exists a sequence of nodes*z*_{1}, … ,*z*_{s}*V*(*G*), with*z*_{1}=*u*and*z*=_{s}*v*, such that (*z*_{1},*z*_{2}), (*z*_{2},*z*_{3}), … , (*z*_{s−1},*z*)_{s}*E*(*G*)). - The descendant nodes $\text{d}{\text{e}}_{w}^{\overline{G}}=\left\{v\in V:{\text{T}}_{wv}^{\overline{G}}=1\right\}$ and the non-descendant nodes $\text{n}{\text{d}}_{w}^{\overline{G}}=\left\{v\in V:{\text{T}}_{wv}^{\overline{G}}=0\right\}$ of
*w**W*in the reduced graph $\overline{G}$. Note that $w\in \text{d}{\text{e}}_{w}^{\overline{G}}\phantom{\rule{0.2em}{0ex}}\text{and}\phantom{\rule{0.2em}{0ex}}w\notin \text{n}{\text{d}}_{w}^{\overline{G}}$ by definition. - The nodes $\text{n}{\text{d}}^{\overline{G}}={\cap}_{w\in W}\text{n}{\text{d}}_{w}^{\overline{G}}$ that are not descendants in the reduced graph $\overline{G}$ of any node in
*W*. - The nodes $\text{d}{\text{e}}_{w}^{\overline{G},H}={\cup}_{x\in \text{p}{\text{a}}_{w}^{H}}\text{d}{\text{e}}_{x}^{\overline{G}}$ that are descendants in the reduced graph $\overline{G}$ of any node $x\in \text{p}{\text{a}}_{w}^{H},$ for a given node
*w**W*. Each node $x\in \text{p}{\text{a}}_{w}^{H}$ is a parent node of the node*w*in the graph $H=\langle {\text{pa}}_{{w}_{1}}^{H},\dots ,{\text{pa}}_{{w}_{q}}^{H}\rangle .$ - The nodes ${\text{de}}_{-w}^{\overline{G},H}={\cup}_{x\in W\backslash {\text{pa}}_{w}^{H}}{\text{de}}_{x}^{\overline{G}}$ that are descenda nts in the reduced graph $\overline{G}$ of any node $x\in W\backslash {\text{pa}}_{w}^{H},$ for a given node
*w**W*. Each node $x\in W\backslash {\text{pa}}_{w}^{H}$ is not a parent node of*w*in the graph $H=\langle \text{p}{\text{a}}_{{w}_{1}}^{H},\dots ,\text{p}{\text{a}}_{{w}_{q}}^{H}\rangle .$

An illustrative example of the relevant graphs and sets, with *W* = {3, 4, 7} shown in red. The edges into *W* are removed from the original graph *G*, to form $\overline{G}$. Here is the set of all DAGs on the nodes {3, 4, 7} and thus *η* = *||* **...**

The relevant sets and requirements of conditions (A) and (B) for the illustrative graph *G* shown in Figure 3, with *H* as shown. Recall ${\text{nd}}^{\overline{G}}=\left\{1,2\right\}$. Condition (B) depends on ${R}_{7,4}^{\overline{G},H}=\left\{4\right\},$ but it imposes no restriction when *w* {3, **...**

Using this notation, we define $P{\text{a}}_{w}^{G,H},$ for given *G* $\mathcal{G}$, *W* *V*, *w* *W* and $H=\langle {\text{pa}}_{{w}_{1}}^{H},\dots ,{\text{pa}}_{{w}_{q}}^{H}\rangle \in \mathscr{H},$ as the set of parent sets ${\text{pa}}_{w}^{{G}^{\prime}}$ that satisfy the following conditions.

$(\text{A})\phantom{\rule{0.5em}{0ex}}{\text{pa}}_{w}^{{G}^{\prime}}\subseteq {Q}_{w}^{\overline{G},H}=\left({\text{nd}}^{\overline{G}}\cup {\text{de}}_{w}^{\overline{G},H}\right)\backslash {\text{de}}_{-w}^{\overline{G},H}$

$(\text{B})\phantom{\rule{0.5em}{0ex}}{\text{pa}}_{w}^{{G}^{\prime}}\cap {R}_{w,x}^{\overline{G},H}\ne \varnothing $ for all nodes $x\in {\text{pa}}_{w}^{H},$ with ${R}_{w,x}^{\overline{G},H}={\text{de}}_{x}^{\overline{G}}\backslash {\text{de}}_{-w}^{\overline{G},H}.$

Note that (B) depends on ${\text{de}}_{x}^{\overline{G}}$ not ${\text{de}}_{x}^{\overline{G},H}.$ We define $P{\text{a}}_{W}^{G,H}$ as the set of tuples formed by the Cartesian product of the sets $P{\text{a}}_{W}^{G,H}$ of parent sets for *w* *W* ; in other words ${\text{pa}}_{W}^{{G}^{\prime}}=({\text{pa}}_{{w}_{1}}^{{G}^{\prime}},\dots ,{\text{pa}}_{{w}_{q}}^{{G}^{\prime}})\in P{\text{a}}_{W}^{G,H}$ if and only if ${\text{pa}}_{{w}_{1}}^{{G}^{\prime}}\in P{\text{a}}_{{w}_{1}}^{G,H},\dots ,{\text{pa}}_{{w}_{q}}^{{G}^{\prime}}\in P{\text{a}}_{{w}_{q}}^{G,H}.$

**Lemma 1**
$\{P{\text{a}}_{W}^{G,{H}_{1}},\dots ,P{\text{a}}_{W}^{G,{H}_{\eta}}\}$
*is a partition of*
$P{\text{a}}_{W}^{G}.$

**Proof** See Appendix B.

Condition (A) ensures all graphs $G\prime =\langle {\text{pa}}_{W}^{G\prime},{\text{pa}}_{-W}^{G}\rangle ,$ with parents ${\text{pa}}_{W}^{G\prime}\in P{\text{a}}_{W}^{G,H},$ are acyclic. It requires that all parents of *w* *W* are either not descendants in the reduced graph $\overline{G}$ of any node in *W*, or a node whose ancestors in the reduced graph $\overline{G}$ are all parents in the graph *H* of node *w*. In particular, no descendant of *w* is added as an ancestor of *w*.

Condition (B) ensures each DAG *G*′ is a member of $P{\text{a}}_{W}^{G,H}$ for only a single *H* for any given *G* $\mathcal{G}$ and thus that $\left\{P{\text{a}}_{W}^{G,{H}_{1}},\dots ,P{\text{a}}_{W}^{G,{H}_{\eta}}\right\}$ is a partition of the set of permissible tuples of parent sets. The condition checks that each edge in the graph *H* is ‘used’ i.e. that ${\text{pa}}_{W}^{G\prime}=\left({\text{pa}}_{{w}_{1}}^{G\prime},\dots ,{\text{pa}}_{{w}_{q}}^{G\prime}\right)\in P{\text{a}}_{W}^{G,H}$ would not be allowed under any *H*′ ≠ *H*, *H*′ . Specifically, it ensures that each parent set ${\text{pa}}_{w}^{G\prime}$ contains at least one descendant node *v* *V* of each parent ${\text{pa}}_{w}^{H}$ of node *w* in the graph *H*, and that *v* is not a descendant in $\overline{G}$ of any node *x* *W* that is not a parent in *H* of *w* *W*. To see why this condition is required, consider a graph *G*′ in which all nodes in *W* have no parents. Then ${\text{pa}}_{w}^{G\prime}=\varnothing $ for all nodes *w* *W* and so condition (A) holds for all *H* . Thus if condition (B) is disregarded, ${\text{pa}}_{W}^{G\prime}\in P{\text{a}}_{W}^{G,H}$ for all *H* , implying that $\left\{P{\text{a}}_{W}^{G,{H}_{1}},\dots ,P{\text{a}}_{W}^{G,{H}_{\eta}}\right\}$ is not a partition of $P{\text{a}}_{W}^{G}.$

Parent sets in each set $P{\text{a}}_{w}^{G,H}$ can be identified easily using simple set operations, and consequently the partition components $P{\text{a}}_{w}^{G,H}$ can be easily identified via Cartesian products of these sets. Let $P{\text{a}}_{w}^{\text{all}}=\mathcal{P}\left(V\backslash \left\{w\right\}\right)$ denote the set of all possible parent sets of node *w*, subject to maximum in-degree *κ*, and $P{\text{a}}_{w}^{\text{all}}(v)=P{\text{a}}_{w}^{\text{all}}\backslash \mathcal{P}\left(V\backslash \left\{w,v\right\}\right)$ denote the subset all $P{\text{a}}_{w}^{\text{all}}$ containing only parent sets that contain node *v* *V*. Parent sets in $P{\text{a}}_{w}^{G,H}$ must (A) not include any nodes *not* in ${Q}_{w}^{\overline{G},H},$ and (B) include at least one node in ${R}_{w,x}^{\overline{G},H}$ for each $x\in {\text{pa}}_{w}^{H},$ and thus

$$P{\text{a}}_{w}^{G,H}=\{\begin{array}{cc}P{\text{a}}_{w}^{(B)}\backslash P{\text{a}}_{w}^{(A)}& {\text{if}\phantom{\rule{0.3em}{0ex}}\text{pa}}_{w}^{H}\ne \varnothing \\ P{\text{a}}_{w}^{\text{all}}\backslash P{\text{a}}_{w}^{(A)}& {\text{if}\phantom{\rule{0.3em}{0ex}}\text{pa}}_{w}^{H}=\varnothing ,\end{array}$$

where $P{\text{a}}_{w}^{(A)}={\displaystyle {\cup}_{v\notin {Q}_{w}^{\overline{G},H}}P{\text{a}}_{w}^{\text{all}}\left(v\right)}\phantom{\rule{0.3em}{0ex}}\text{and}\phantom{\rule{0.3em}{0ex}}P{\text{a}}_{w}^{(B)}={\displaystyle {\cap}_{x\in {\text{pa}}_{w}^{H}}{\displaystyle {\cup}_{r\in {R}_{w,x}^{\overline{G},H}}P{\text{a}}_{w}^{\text{all}}\left(r\right)}}.$

We fix *q* = *|W|* as a small constant for all *p* and set a maximum in-degree *κ*. This means permissible tuples of parents sets can be identified in $\mathcal{O}$(*p*^{κ+1}) time by storing the ‘lookup tables’ $P{\text{a}}_{w}^{\text{all}}\left(v\right)$ as bit maps, assuming that $\mathcal{O}$(1) querying of the descendants and non-descendants is available.

This can be achieved using a dynamic algorithm that provides access to the transitive closure. Giudici and Castelo (2003) previously used a similar approach. For a graph with transitive closure T^{G} with edges *E*^{T}, the descendants and non-descendants of node *u* are {*v* : (*u*, *v*) *E*^{T}} and *{v* : (*u*, *v*) *E*^{T}*}* respectively. Thus, the descendants and non-descendants of a node can be identified in $\mathcal{O}$(1) time when the transitive closure is known. The transitive closure for an arbitrary directed graph with *p* vertices can be determined in $\mathcal{O}$(*p ^{ω}*) time (Munro, 1971), where

The algorithm maintains a *p* × *p* path count matrix *C ^{G}* whose elements ${C}_{uv}^{G}$ are the number of distinct paths from node

We draw a new graph $G\prime =\langle {\text{pa}}_{W}^{G\prime},{\text{pa}}_{-W}^{G}\rangle $ starting from a graph *G* using a two-stage method: we sample first a component $P{\text{a}}_{W}^{G,H\prime}$ of the partition of permissible tuples of parent sets and then a tuple ${\text{pa}}_{W}^{G\prime}$ of parent sets from the selected partition component.

In the first stage, $P{\text{a}}_{W}^{G,H\prime}$ is drawn from the conditional distribution, given the tuple of parent sets of nodes not in *W*.

$$\begin{array}{l}P(P{\text{a}}_{W}^{G,H\prime}|{\text{pa}}_{-W}^{G},\mathbf{\text{X}})=\frac{P\left(P{\text{a}}_{W}^{G,H\prime},{\text{pa}}_{-W}^{G}|\mathbf{\text{X}}\right)}{P\left({\text{pa}}_{-W}^{G}|\mathbf{\text{X}}\right)}\\ \phantom{\rule{8.2em}{0ex}}=\frac{{\displaystyle {\sum}_{{\text{pa}}_{W}^{G\prime}\in P{\text{a}}_{W}^{G,H\prime}}{\displaystyle {\prod}_{w\in W}p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}^{G\prime}}\right){\pi}_{w}\left({\text{pa}}_{w}^{G\prime}\right)}}}{{\displaystyle {\sum}_{H\u2033\in \mathscr{H}}{\displaystyle {\sum}_{{\text{pa}}_{W}^{G\u2033}\in {\text{Pa}}_{W}^{G,H\u2033}}{\displaystyle {\prod}_{w\in W}p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}^{G\u2033}}\right){\pi}_{w}\left({\text{pa}}_{w}^{G\u2033}\right)}}}}\\ \phantom{\rule{8.2em}{0ex}}=\frac{{\displaystyle {\prod}_{w\in W}{\displaystyle {\sum}_{{\text{pa}}_{w}^{G\prime}\in P{\text{a}}_{w}^{G,H\prime}}p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}^{G\prime}}\right){\pi}_{w}\left({\text{pa}}_{w}^{G\prime}\right)}}}{{\displaystyle {\sum}_{H\u2033\in \mathscr{H}}{\displaystyle {\prod}_{w\in W}{\displaystyle {\sum}_{{\text{pa}}_{w}^{G\u2033}\in P{\text{a}}_{w}^{G,H\u2033}}p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}^{G\u2033}}\right){\pi}_{w}\left({\text{pa}}_{w}^{G\u2033}\right)}}}}\end{array}$$

The final equality follows by an interchange of sum and product that is proved in Lemma 2. This makes evaluation more efficient by allowing the sums to be evaluated separately for each node. Friedman and Koller (2003) used a similar interchange.

**Lemma 2**
*The following identity holds for any H* , *W* *V*
*and*
*G* $\mathrm{\mathcal{G}.}$

$$\sum _{{\text{pa}}_{W}\in P{\text{a}}_{W}^{G,H}}{\displaystyle \prod _{w\in W}p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}}\right){\pi}_{w}\left({\text{pa}}_{w}\right)}}={\displaystyle \prod _{w\in W}{\displaystyle \sum _{{\text{pa}}_{w}\in P{\text{a}}_{w}^{G,H}}p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}}\right){\pi}_{w}\left({\text{pa}}_{w}\right)}$$

**Proof** See Appendix C.

In the second stage, we sample new parents ${\text{pa}}_{W}^{{G}^{\prime}}$ from the selected partition component, and form the new graph ${G}^{\prime}=\langle {\text{pa}}_{W}^{{G}^{\prime}},{\text{pa}}_{-W}^{G}\rangle .$ The parents of each node *w* *W*, conditional on *H*′, are independent, and so can be sampled separately from the following conditional distribution:

$$P\left({\text{pa}}_{w}^{{G}^{\prime}}|P{\text{a}}_{W}^{G,{H}^{\prime}},{\text{pa}}_{-W}^{G},\mathbf{\text{X}}\right)=\frac{p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}^{{G}^{\prime}}}\right){\pi}_{w}\left({\text{pa}}_{w}^{{G}^{\prime}}\right)}{{\displaystyle {\sum}_{{\text{pa}}_{w}^{{G}^{\u2033}}\in P{\text{a}}_{w}^{G,{H}^{\prime}}}p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}^{{G}^{\u2033}}}\right){\pi}_{w}\left({\text{pa}}_{w}^{{G}^{\u2033}}\right)}}.$$

This step is straightforward because this distribution is simply the posterior distribution of a constrained variable selection with response *w* and $P{\text{a}}_{w}^{G,{H}^{\prime}}$ as the set of possible active sets (i.e. selected covariates).

Algorithm 2 outlines the complete algorithm. The methods in Sections 4.2 enable fast identification of each partition component. Run-time is a function of *p*, maximum in-degree *κ*, and the number *q* of nodes in the block. We choose a small, fixed *q* for all *p*, so the run-time is determined by the evaluation of $P{\text{a}}_{w}^{G,H},$ which is $\mathcal{O}$(*p*^{κ+1}).

Initialise starting point ${G}_{0}=\langle {\text{pa}}_{1}^{{G}_{0}},\dots ,{\text{pa}}_{p}^{{G}_{0}}\rangle $

Compute initial path count matrix *C*^{G0}

**for**
*t* in 1 to *N*
**do**

Sample *q* nodes uniformly at random from *V* , and call this set of nodes *W*

Let *G* = *G*_{t−1}

Form $\overline{G}$ as defined in Section 4.1

Evaluate ${C}^{\overline{G}}$ from *C ^{G}* as described in Section 4.2

**for**
*w* *W*
**do**

Evaluate ${\text{de}}_{w}^{\overline{G}}=\left\{v:\phantom{\rule{0.3em}{0ex}}{C}_{wv}^{\overline{G}}\ge 1\right\}$

Evaluate ${\text{nd}}_{w}^{\overline{G}}=\left\{v:\phantom{\rule{0.3em}{0ex}}{C}_{wv}^{\overline{G}}=0\right\}$

**end for**

**for**
*H* **do**

**for**
*w* *W*
**do**

Evaluate $P{\text{a}}_{w}^{G,H},$ as described in Section 4.2

Let ${K}_{w}^{H}={\displaystyle {\sum}_{{\text{pa}}_{w}^{{G}^{\prime}}\in P{\text{a}}_{w}^{G,H}}p\left({\mathbf{\text{X}}}_{w}|{\mathbf{\text{X}}}_{{\text{pa}}_{w}^{{G}^{\prime}}}\right){\pi}_{w}\left({\text{pa}}_{w}^{{G}^{\prime}}\right)}$

**end for**

Let ${K}^{H}={\displaystyle {\prod}_{w\in W}{K}_{w}^{H}}$

**end for**

Sample $P{\text{a}}_{W}^{G,{H}^{\prime}}$ according to $P\left(P{\text{a}}_{W}^{G,{H}^{\prime}}|{\text{pa}}_{-W}^{G},\mathbf{\text{X}}\right)=\frac{{K}^{{H}^{\prime}}}{{\displaystyle {\sum}_{{H}^{\u2033}\in \mathscr{H}}{K}^{{H}^{\u2033}}}}$

**for**
*w* *W*
**do**

Sample ${\text{pa}}_{W}^{{G}^{\prime}}$ according to $P\left({\text{pa}}_{w}^{{G}^{\prime}}|P{\text{a}}_{W}^{G,{H}^{\prime}},{\text{pa}}_{-W}^{G},\mathbf{\text{X}}\right)$

**end for**

Set ${G}_{t}\leftarrow G=\langle {\text{pa}}_{W}^{{G}^{\prime}},{\text{pa}}_{-W}^{{G}_{t-1}}\rangle $

Update *C*^{Gt}

**end for**

In this section, we empirically assess the performance of the proposed sampler, comparing it with existing samplers as well as frequentist methods for structure learning of DAGs.

We compare the Gibbs sampler with the MC^{3} sampler (Madigan and York, 1995) and the REV sampler (Grzegorczyk and Husmeier, 2008), a variant of MC^{3} that uses a more extensive edge reversal move. We also compare with two frequentist constraint-based methods: the PC-algorithm (Spirtes et al., 2000), and the Xie and Geng (2008) method, that is shown by its authors to outperform the PC-algorithm in some settings.

Tuning parameters for each method were set as follows. For the constraint-based methods, the significance level was *α* = 0.05 by default, but we also show some results for *α* = 0.00005, 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.1. The Gibbs sampler we use is a random-scan sampler, with *q* = 3 (i.e. the parent sets of three nodes are sampled jointly at each iteration). To permit a fair comparison, for MC^{3} we use the same fast online transitive closure algorithm (Section 4.2), and pre-computation and caching of local marginal likelihoods used in the Gibbs sampler (REV uses a similar pre-computation and caching scheme). We constrain all of the samplers to maximum in-degree *κ* = 3 and set the graph prior as *π*(*G*) 1. We use conjugate formulations throughout, specifically multinomial-Dirichlet (Heckerman et al., 1995) for discrete data and Normal with a *g*-prior (with *g* = *n*) (Geiger and Heckerman, 1994; Zellner, 1986) for continuous data.

We consider six examples: a small domain example, where comparison with the exact posterior (Tian and He, 2009) is possible; data simulated from the ALARM network and from randomly-generated networks of varying sparsity; data sets from social science and biology; and a pathological 4-node example designed to highlight a failure case for our method.

In practical use samplers can only be run for a finite number of iterations (depending on available time and computational resources). We set the maximum number of iterations as follows. In total, we drew 10^{6} iterations of REV (retaining only every 10^{th} iteration to reduce storage requirements). Following Grzegorczyk and Husmeier (2008), 85% of proposals within REV were MC^{3} proposals (without MC^{3} proposals, the REV sampler is not irreducible). In our implementation, the computational costs of the Gibbs sampler are an order of magnitude lower than REV’s (accounting for MC^{3} moves), but we nevertheless treated the computational costs as the same for the purposes of comparison and drew 10^{6} iterations of the Gibbs sampler (again retaining every 10^{th} iteration). The computational costs of our implementation of MC^{3} are roughly 1/10 of a Gibbs iteration, and so we performed 10^{7} iterations of MC^{3} (retaining every 100^{th}).

For each sampler, 10 independent runs starting from different initial graphs were performed. We discarded the first 1/4 of samples, but as an additional filter we considered trace plots (log marginal likelihood versus iteration) for each run of each sampler and discarded further samples if they appeared to be from a pre-convergent phase. Specifically, if there was a clear ‘step change’ in the trace for a certain run we discarded all samples in that run drawn prior to the highest ‘step’ being reached. We note that this simple filter cannot be regarded as detecting convergence because the highest ‘step’ may correspond to a region near a local rather than a global mode. We therefore also study convergence in detail via other metrics.

Each sampler gives an estimated posterior distribution over DAGs. From this we obtain a point estimate either as the **maximum a posteriori** (MAP) graph *G*^{MAP}, or as a **thesholded graph**
*G _{τ}* formed by including all directed edges whose marginal posterior edge probabilities are at least

We use the structural Hamming distance (SHD) to quantify differences between DAGs. This is the minimum number of edge insertions and deletions needed to transform one graph into another. Comparisons are made using completed partially directed acyclic graphs (Chickering, 2002a). To show the behaviour of the samplers at shorter MCMC run lengths some metrics are shown against iteration *t*. These are based on the final 3/4 of samples drawn up to *t* (except for log marginal likelihoods), and so may incorporate some samples drawn before convergence.

We assess convergence and stability via the following metrics:

**Trace plots**of log marginal likelihood against iteration*t*.- Between-run agreement between
**posterior edge probabilities**. These are visualized as scatter plots. When edge probabilities agree between a pair of runs, all points in the plot will lie close to the*y*=*x*line. We show two variants: a hexagonally-binned version, to avoid over-plotting (Carr et al., 1987); and a panelled plot, in which each panel shows one pair of runs. We also consider the number of edges with estimated posterior probability greater than 0.9 in one run, and less than 0.1 in another run. Such**‘major discrepancies’**represent serious Monte Carlo artefacts. **Potential scale reduction factor (PSRF)**is a multiple-chain convergence diagnostic, developed by Gelman and Rubin (1992), that compares within- and between-chain means and variances of a suitable summary statistic. PSRF < 1.1 is often taken to indicate that a sampler has run sufficiently long. The natural summary statistics in this context are the posterior edge probabilities. However, we find that the resulting PSRF is not a reliable indicator of convergence. In simulations (not shown) using a 20 node network for which we could calculate true posterior edge probabilities (Tian and He, 2009), we found little association between single-edge PSRF and the (absolute) error in posterior edge probability. In particular, the diagnostic appeared to be very conservative for edges with true posterior probability near 0 or 1. We therefore suggest that a reasonable use of PSRF in this context may be to assess relative convergence between different samplers, rather than as an absolute indicator of convergence. Below, we compare samplers in terms of the proportion of edges with PSRF < 1.1 (a larger proportion suggests better mixing). To calculate PSRF we consider the 10 runs of each sampler as a collection of 5 pairs of runs and calculate PSRF separately for each edge using the final three quarters of the samples drawn up to that point.- For the real data, we assess
**stability under resampling**(“shaking the data”) by comparing estimates across bootstrap samples.

For experiments in which the true data-generating graph *G** is known we use the following metrics to assess accuracy:

**Structural Hamming distance**(SHD) between*G** and the estimated graphs.**Receiver-operating characteristic (ROC) curves**. These show agreement between*G** and an estimate*G*by plotting true positive against false positive rates parameterized by threshold_{τ}*τ*. We consider also the**area under the ROC curve**(AUROC), focusing in particular on the small false positive rate region of the curve that is often of interest in applications.

Finally, when the posterior edge probabilities can be computed exactly (Tian and He, 2009) (i.e. in small *p* settings) we also consider the **maximum and average absolute error in posterior edge probability**, calculated across the set of all possible edges.

We note that while SHD and ROC scores are useful in assessing performance, they are not convergence measures, as they do not assess accuracy of the posterior distribution *per se* (to see why, consider a degenerate sampler that samples only *G**, giving perfect scores on SHD and ROC, but an incorrect posterior distribution).

We applied the methods to the Zoo data set (Newman et al., 1998) that records *p* = 17 (discrete) characteristics of *n* = 101 animals. The maximum log marginal likelihood found by the Gibbs and REV samplers was consistent across runs, but less so for the MC^{3} sampler, and the Gibbs sampler reached a plateau of high probability after the fewest iterations (Figure A13a). REV required about ten times as many iterations as Gibbs to achieve the same proportion of edges with PSRF < 1.1 while MC^{3} needed about 100 times as many (Figure A14a). The estimated edge probabilities given by the Gibbs sampler were stable between runs (Figure A15a). The results from the REV sampler were also stable, but MC^{3} less so, with major disparities between some runs. Figure 4 shows the error in posterior edge probabilities as a function of iterations. Convergence was quickest for the Gibbs sampler, followed by REV and then MC^{3}. All runs of the Gibbs sampler reached a point at which the maximum absolute error was 0.05 (after 67,000 iterations on average). In contrast, only 5/10 runs of REV and 6/10 runs of MC^{3} reached the same level of accuracy in their complete run. Similarly, the Gibbs sampler achieved an average absolute error of less than 0.01 in the fewest iterations.

The ALARM network (Beinlich et al., 1989) consists of 37 discrete nodes and 46 edges and has been widely used in studying structure learning (e.g. Friedman and Koller, 2003; Grzegorczyk and Husmeier, 2008). We simulated data sets from ALARM with sample sizes *n* = 100, 500, 1000, 2500, 5000. Figure 5 shows trace plots for the *n* = 1000 case as an illustrative example. The maximum log marginal likelihood found by the Gibbs sampler across runs was −10,608.20 ± 0.8 (mean ± standard deviation), whereas for REV it was −10,627.2 ± 40.4 and for MC^{3} −11,176.9 ± 273.7. The highest scoring graph discovered by any of the 10 runs of the MC^{3} sampler had log marginal likelihood −10,856.6, far below the Gibbs maximum, suggesting non-convergence in all 10 runs. REV appeared to reach convergence in all but two runs (runs 9 and 10). The Gibbs sampler consistently (and rapidly) reached a high scoring plateau and appeared to have converged in all 10 runs. PSRF results were in line with these findings (Figure A14b); more than 10 times as many iterations of REV were needed to give the same proportion of edges with PSRF < 1.1 as the Gibbs sampler.

Log marginal likelihood of the graphs visited by the three MCMC samplers in 10 independent runs, initialised at disparate initial conditions, on the ALARM data, with data sample size *n* = 1000. Iteration number is displayed on a log_{10} scale. The dashed **...**

Figure 6 compares pairs of runs for *n* = 1000; this pair of runs is typical of all 10 runs (except for runs 9 and 10 of REV which disagree considerably; all runs shown in Figure A15b). There were no major between-run discrepancies (as defined in Section 5.2) for the Gibbs sampler at any sample size. The mean number of major discrepancies (across pairs of runs) increased from 0 (*n* = 100) to 8 and 91 for MC^{3} and REV respectively (when *n* = 5000).

Convergence diagnostic plots for each MCMC sampler for the ALARM data, with *n* = 1000. The posterior edge probabilities given by two independent MCMC runs are plotted against each other, binned into hexagonal areas to avoid over-plotting. When the edge **...**

Figure 7 shows ROC curves for false positive rates < 0.05. The Gibbs sampler performs better and with less variability than the other methods. Sample size is influential. For small *n* the Bayesian methods outperform the constraint-based methods. However, counter to the increase in statistical information, REV and MC^{3} perform less well with larger *n*: e.g. at *n* = 100 for a false positive rate (FPR) of 0, the Gibbs sampler found 21.9 ± 1.9 (mean ± standard deviation) true edges; REV 20.1 ± 1.4; and MC^{3} 21.7 ± 2.1. But at *n* = 5000, for the same FPR, Gibbs found 43.0 ± 0.0 true edges; REV 13.2 ± 21.3; and MC^{3} did not find any true edges (the true graph has 46 edges). The constraint-based methods performed well for large sample sizes, as anticipated by the asymptotic consistency of the PC-algorithm (Kalisch and Bühlmann, 2007). The Xie-Geng method performed particularly well for *n* = 5000. Figure 8 shows AUROC as a function of iterations; the Gibbs sampler performed best after all run lengths and sample sizes considered. These results are supported by SHDs between point estimates and the true graph (Table A2).

ALARM data, receiver operating characteristic (ROC) curves for each of 10 independent MCMC runs for each MCMC sampler, and point estimates for the constraint-based methods. Note that the horizontal axis shows only false positives rates < 0.05, **...**

Sparsity can be used to statistical and computational advantage but in practice it may be hard to know what level of sparsity is reasonable for a given application. We therefore sought to investigate the effect of varying network sparsity, including scenarios where the data-generating graph can violate the in-degree constraint we impose. We simulated data following a procedure described in Kalisch and Bühlmann (2007) that we outline below. We first generated a DAG *G* via its adjacency matrix *A ^{G}*, by drawing entries as ${A}_{uv}^{G}\sim \text{Bernoulli(}\rho \text{),}$ where

Figure 9 shows boxplots over SHDs. We see that accuracy decreases with increasing density, echoing results in Kalisch and Bühlmann (2007) for the PC algorithm. As throughout, all the samplers had an in-degree constraint (maximum in-degree *κ* = 3), while the frequentist methods did not. Nevertheless, the performance of the Bayesian methods did not appear to deteriorate any more rapidly than the frequentist methods. At all *ρ*’s, the median probability graph *G*^{med} outperformed the frequentist methods which in turn outperformed the MAP graph. In this context, we draw attention to a difference between the median probability and MAP graphs: the former, although obtained by averaging over DAGs satisfying the in-degree constraint, may itself have in-degree greater than *κ*, while the latter is necessarily subject to the constraint.

The publicly available Behavioral Risk Factor Surveillance System Survey (BRFSS) (Centers for Disease Control and Prevention, 2008) is a household-level random-digit telephone survey, collected by the U.S. National Center for Chronic Disease Prevention and Health, that has been conducted throughout the United States since 1984. We considered (discrete) responses to *p* = 24 questions (see Appendix D), spanning most of the topics covered in the survey. We considered the responses from New York in the 2008 survey, removing samples for respondents who refused or were unsure of their response, or whose response was missing, to any of the 24 questions. The resulting sample size was *n* = 4, 197. The median probability graph *G*^{med} estimated by the Gibbs sampler is shown in Figure A17a.

Figure 10 shows between-run agreement. The Gibbs runs showed better agreement than the REV runs and the MC^{3} runs disagreed considerably (these pairs of runs were typical; all runs shown in Figure A18a). Indeed there were no major between-run discrepancies (in the sense of Section 5.2) for the Gibbs sampler, whereas there were on average 2.4 major discrepancies for REV and 10.9 for MC^{3}. The Gibbs sampler also had the highest proportion of edges with PSRF < 1.1 (Figure A14c).

Convergence diagnostic plots for each MCMC sampler for the survey data. The posterior edge probabilities given by two independent MCMC runs are plotted against each other, binned into hexagonal areas to avoid over-plotting. When the two runs give the **...**

The maximum log marginal likelihoods found by each of the three samplers differed considerably as did the number of iterations needed to reach a plateau (Figure A13b). The Gibbs sampler typically reached a plateau after around 500 samples (although in one run 10^{3} samples were needed). REV took longer to reach a (usually lower) plateau. MC^{3} appeared to become stuck in a region with even lower log marginal likelihood.

To investigate stability under resampling of the data we applied the methods to 10 bootstrap resamples of the data set. The Gibbs and REV samplers were more stable than the frequentist methods as well as than MC^{3}. For example, using ${G}_{\tau}^{\text{PC}}$ as a point estimate for for the Bayesian methods (i.e. thresholding to give the same number of edges as PC), the mean SHD between results from pairs of bootstrap data sets was 21.7 (Gibbs), 22.3 (REV), 33.4 (MC^{3}) and 30.8 (PC-algorithm). Results for *G*^{MAP} and ${G}_{\tau}^{\text{Xie}}$ are shown in Figure A19a.

We used single-cell molecular data from Bendall et al. (2011) with *p* = 34 continuous variables (see Appendix D) and *n* = 21, 691. The median probability graph *G*^{med} estimated by the Gibbs sampler is shown in Figure A17b.

Figure 11 shows between-run agreement for the samplers for a typical pair of runs; the Gibbs sampler shows better agreement than REV or MC^{3}. All but one of the 10 runs of the Gibbs sampler were consistent with each other (Figure A18b) and the one inconsistent run nonetheless showed better agreement with the other Gibbs runs than any pair of runs of the other samplers. The Gibbs sampler had no major discrepancies (as defined in Section 5.2) between any pairs of runs, while there were on average 26 and 87 major discrepancies for REV and MC^{3} respectively. Smaller discrepancies followed a similar pattern: the average number of edges differing by 0.1 or more in posterior probability between runs were 6.4 (Gibbs sampler), 56.6 (REV), and 151.8 (MC^{3}). The Gibbs sampler also had the highest proportion of edges with PSRF < 1.1 after any run length (Figure A14d). Trace plots are shown in Figure A13c. The Gibbs sampler found a region of higher log marginal likelihood than the other samplers and did so consistently. Indeed, the maximum log marginal likelihood achieved across all REV runs was exceeded in every Gibbs run and after only around 5000 iterations. Finally, we considered stability under resampling of the data (Figure A19b), including also the frequentist point estimators for comparison. The mean SHD between PC estimates from pairs of bootstrap samples was 214.6; the corresponding mean SHDs for the samplers (using ${G}_{\tau}^{\text{PC}}$ point estimates) were 128.8 (Gibbs), 140.1 (REV) and 225.2 (MC^{3}). We note that the SHDs are particularly high here because the PC-algorithm estimates a network with a high density (recall that under ${G}_{\tau}^{\text{PC}}$ all methods choose the same number of edges as PC).

Our final example demonstrates the potential sensitivity of the Gibbs sampler to choice of *q* = *|W|* (the number of nodes whose parent sets are sampled together in a single iteration). The example was constructed to highlight a situation in which the sampler can become stuck in a local mode unless *q* is large enough, potentially leading to extremely slow convergence.

We used as data 10^{5} simulated samples from a 4-node network in which the parents of node 2 were nodes 1 and 4; the parents of node 3 were nodes 1 and 2; and nodes 1 and 4 had no parents. Each node was Bernoulli distributed. The probability of success was 0.6 for nodes 1 and 4, while nodes 2 and 3 were noisy XORs with probability of success 0.9 when either parent (but not both) was ‘true’ and 0.1 otherwise.

For the purposes of demonstration, we set *q* = 2. As shown in Figure 12, after 10^{6} iterations the maximum error in edge probability for the Gibbs sampler was 0.016 *±* 0.015 (mean *±* standard deviation across 10 runs). Given that there are only 543 DAGs with *p* = 4 nodes, this magnitude of error is unexpectedly large. The slow convergence is due to the concentration of posterior mass on two graphs that the sampler cannot easily move between. The MAP graph (posterior probability 0.61) is the graph in which the parents of node 2 are nodes 1 and 3; the parents of node 4 are nodes 1 and 2; and nodes 1 and 3 have no parents. The data-generating graph has posterior probability 0.38. The graphs differ in the parents of 3 nodes, and so with *q* = 2 the sampler cannot move between them in a single step. At the same time, the large sample size leads to all other graphs having low posterior probability (< 0.002), making multi-step transitions between the two graphs unlikely. Thus, the sampler becomes stuck on one of the two graphs. In this example, setting *q* = 3 improved convergence (maximum error 0.0011 *±* 0.0009 after 10^{6} iterations), but in real-world examples it is difficult to rule out the possibility that analogous issues may arise.

We introduced a Gibbs sampler for structure learning of DAGs that can converge more easily than existing samplers due to its ability to efficiently make large moves in DAG space. We showed that it provides often substantial gains in accuracy and stability in comparison with existing (Bayesian and frequentist) methods in a range of settings.

The formulation of the sampler develops and exploits the connection between variable selection and structure learning. This connection has been widely studied for undirected graphs (e.g. Meinshausen and Bühlmann, 2006), but for DAGs is complicated by the acyclicity requirement. Our approach accounts for acyclicity but further work in this area may ease the adaptation of results and methods from Bayesian variable selection to the case of DAGs. In the proposed sampler, existing variable selection methods could be of direct utility in sampling from the conditional distribution $P(\text{p}{\text{a}}_{W}|\text{p}{\text{a}}_{-W}^{G},\mathbf{\text{X}}).$ When this sampling step is difficult, a Metropolis-within-Gibbs approach (i.e. substituting a Metropolis step in place of the Gibbs step) could be considered. With *q* = *|W|* = 1, the conditional $P(\text{p}{\text{a}}_{W}|\text{p}{\text{a}}_{-W}^{G},\mathbf{\text{X}})$ is identical to the posterior of the corresponding variable selection problem. Then, a Metropolis-within-Gibbs move could directly exploit existing variable selection methods.

In common with most of the structure learning literature, we used an in-degree constraint. This gives a smaller DAG space and in addition controls the number of parameters needed to specify conditional distributions. However, it is a strong constraint and a natural concern is whether it excludes higher in-degree models that could be appropriate for the data. This possibility cannot be ruled out, but the use of model averaging and marginal posterior summaries may ameliorate the concern to some extent, because even when the true number of parents of a node exceeds the restriction, important candidate parents may still have high edge scores (and appear in a summary graph). For example, suppose the true in-degree of a node is 4, but at most 3 incoming edges are permitted. In this case, although no model including all 4 parents can be considered, provided the signal can be detected considering only 3 nodes at a time, the posterior probability of all 4 edges may nonetheless be high.

In our empirical examples the REV sampler of Grzegorczyk and Husmeier (2008) showed impressive gains over MC^{3}. The Gibbs sampler generally seemed to outperform both. The use of conditional distributions in REV is a point of similarity with the Gibbs sampler proposed here. The performance gains of Gibbs could be explained by two key differences from REV. First, REV does not use the natural conditional distribution and requires an accept-reject step. Second, REV must include at least some MC^{3} proposals (otherwise the sampler is not irreducible), and these steps are not tailored to the shape of the posterior distribution (Grzegorczyk and Husmeier, 2008, make REV proposals with probability 1/15 and so most steps are based on MC^{3} proposals).

If there is strong correlation between the parent sets of more than *q* = *|W|* nodes, the Gibbs sampler may not mix well. In this situation, constraint-based methods may be useful. Alternatively, *q* could be chosen at each step according to some distribution, so that a mixture of different block sizes is used. This would in particular allow larger blocks to be used without increasing the computational demands of the algorithm excessively. In this paper we fixed *q* = 3, and found this simple choice gave a well-behaved and effective sampler. But there is a trade-off: increasing *q* increases the time taken to evaluate $P(\text{p}{\text{a}}_{W}|\text{p}{\text{a}}_{-W}^{G},\mathbf{\text{X}})$, but also increases move size, with the potential to improve convergence.

Practical use of the Gibbs sampler in the form described here requires exact sampling from the conditional $P(\text{p}{\text{a}}_{W}|\text{p}{\text{a}}_{-W}^{G},\mathbf{\text{X}})$ and there are situations related to this requirement in which other methods may be more suitable. First, when an appropriate maximum in-degree cannot be set, MC^{3} or a variant could be more appropriate (although convergence could be very slow). Alternatively, search procedures such as GES (Chickering, 2002b) or HCMC (Castelo and Kočka, 2003) could be used. Second, exact sampling from the conditional distribution is challenging in settings with thousands of nodes. In this case, efficient constraint-based methods (such as the PC-algorithm) may be a better choice, particularly in the large sample setting. As noted by a referee, an interesting area for future research would be combining the Gibbs sampler with some aspects of other methods—such as PC—that are relatively well suited to the truly high-dimensional setting.

The authors are grateful to Karen Sachs for providing the single-cell molecular data; to Marco Grzegorczyk and Dirk Husmeier for their inspiring work in this area and for providing their implementation of the REV sampler; to the referees for numerous suggestions that improved the article; and to Frances Griffiths, David Lunn and Paul Newcombe for helpful discussions. Part of this work was carried out while the authors were at the University of Warwick—whose excellent scientific environment the authors warmly acknowledge—and supported by the Economic and Social Research Council (ESRC) and Engineering and Physical Sciences Research Council (EPSRC).

Robert J. B. Goudie, Medical Research Council Biostatistics Unit Cambridge CB2 0SR, UK.

Sach Mukherjee, German Centre for Neurodegenerative Diseases (DZNE) Bonn 53175, Germany.

- Barbieri MM, Berger JO. Optimal predictive model selection. Annals of Statistics. 2004;32(3):870–897.
- Beinlich IA, Suermondt HJ, Chavez RM, Cooper GF. Second European Conference on Artificial Intelligence in Medicine. Springer-Verlag; Berlin: 1989. The ALARM monitoring system: A case study with two probabilistic inference techniques for belief networks; pp. 247–256.
- Bendall SC, Simonds EF, Qiu P, Amir ED, Krutzik PO, Finck R, Bruggner RV, Melamed R, Trejo A, Ornatsky OI, Balderas RS, et al. Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum. Science. 2011;332(6030):687–696. [PMC free article] [PubMed]
- Besag JE. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B (Methodological) 1974;36(2):192–236.
- Besag JE. Discussion of “Markov chains for exploring posterior distributions” Annals of Statistics. 1994;22(4):1734–1741.
- Carr DB, Littlefield RJ, Nicholson WL, Littlefield JS. Scatterplot matrix techniques for large N. Journal of the American Statistical Association. 1987;82(398):424–436.
- Castelo R, Kočka T. On inclusion-driven learning of Bayesian networks. Journal of Machine Learning Research. 2003;4:527–574.
- Centers for Disease Control and Prevention. Behavioral Risk Factor Surveillance System Survey Data. U.S. Department of Health and Human Services; Atlanta, Georgia: 2008.
- Chickering DM. Learning equivalence classes of Bayesian-network structures. Journal of Machine Learning Research. 2002a;2:445–498.
- Chickering DM. Optimal structure identification with greedy search. Journal of Machine Learning Research. 2002b;3:507–554.
- Colombo D, Maathuis MH. Order-independent constraint-based causal structure learning. Journal of Machine Learning Research. 2014;15:3921–3962.
- Coppersmith D, Winograd S. Matrix multiplication via arithmetic progressions. Journal of Symbolic Computation. 1990;9(3):251–280.
- Demetrescu C, Italiano GF. Trade-offs for fully dynamic transitive closure on DAGs: breaking through the
*O*(*n*^{2}) barrier. Journal of the ACM. 2005;52(2):147–156. - Demetrescu C, Eppstein D, Galil Z, Italiano GF. Dynamic graph algorithms. In: Atallah MJ, Blanton M, editors. Algorithms and Theory of Computation Handbook: General Concepts and Techniques. Chapman and Hall/CRC; Boca Raton, FL: 2010. pp. 9.1–9.28.
- Eaton D, Murphy K. Bayesian structure learning using dynamic programming and MCMC. Proceedings of the Twenty-Third Conference on Uncertainty in Artificial Intelligence; Corvallis, Oregon: AUAI Press; 2007. pp. 101–108.
- Ellis B, Wong WH. Learning causal Bayesian network structures from experimental data. Journal of the American Statistical Association. 2008;103(482):778–789.
- Friedman N, Koller D. Being Bayesian about network structure. A Bayesian approach to structure discovery in Bayesian networks. Machine Learning. 2003;50(1–2):95–125.
- Geiger D, Heckerman D. Learning Gaussian networks. Proceedings of the Tenth Conference on Uncertainty in Artificial Intelligence; San Francisco, CA: Morgan Kaufmann; 1994. pp. 235–240.
- Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Statistical Science. 1992;7(4):457–472.
- Giudici P, Castelo R. Improving Markov chain Monte Carlo model search for data mining. Machine Learning. 2003;50(1–2):127–158.
- Grzegorczyk M, Husmeier D. Improving the structure MCMC sampler for Bayesian networks by introducing a new edge reversal move. Machine Learning. 2008;71(2–3):265–305.
- Heckerman D, Geiger D, Chickering DM. Learning Bayesian networks: the combination of knowledge and statistical data. Machine Learning. 1995;20(3):197–243.
- Hobert JP, Robert CP, Goutis C. Connectedness conditions for the convergence of the Gibbs sampler. Statistics & Probability Letters. 1997;33(3):235–240.
- Kalisch M, Bühlmann P. Estimating high-dimensional directed acyclic graphs with the PC-Algorithm. Journal of Machine Learning Research. 2007;8:613–636.
- King V, Sagert G. A fully dynamic algorithm for maintaining the transitive closure. Journal of Computer and System Sciences. 2002;65(1):150–167.
- Koivisto M, Sood K. Exact Bayesian structure discovery in Bayesian networks. Journal of Machine Learning Research. 2004;5:549–573.
- Koller D, Friedman N. Probabilistic Graphical Models: Principles and Techniques. MIT Press; Cambridge, MA: 2009.
- Korb KB, Nicholson AE. Bayesian Artificial Intelligence. Chapman & Hall/CRC Press; Boca Raton, FL: 2011.
- Madigan D, York JC. Bayesian graphical models for discrete data. International Statistical Review / Revue Internationale de Statistique. 1995;63(2):215–232.
- Meinshausen N, Bühlmann P. High-dimensional graphs and variable selection with the Lasso. Annals of Statistics. 2006;34(3):1436–1462.
- Munro I. Efficient determination of the transitive closure of a directed graph. Information Processing Letters. 1971;1(2):56–58.
- Newman D, Hettich S, Blake C, Merz C. UCI Repository of Machine Learning Databases. University of California, Department of Information and Computer Science; Irvine, CA: 1998. URL http://www.ics.uci.edu/~mlearn/MLRepository.html.
- Nikolova O, Zola J, Aluru S. Parallel globally optimal structure learning of Bayesian networks. Journal of Parallel and Distributed Computing. 2013;73(8):1039–1048.
- Parviainen P, Koivisto M. Exact structure discovery in Bayesian networks with less space. Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence; Corvallis, Oregon: AUAI Press; 2009. pp. 436–443.
- Parviainen P, Koivisto M. Finding optimal Bayesian networks using precedence constraints. Journal of Machine Learning Research. 2013;14:1387–1415.
- Roberts GO, Sahu SK. Updating schemes, correlation structure, blocking and parameterization for the Gibbs sampler. Journal of the Royal Statistical Society: Series B (Methodological) 1997;59(2):291–317.
- Spirtes P, Glymour C, Scheines R. Causation, Prediction, and Search. The MIT Press; Cambridge, MA: 2000.
- Tamada Y, Imoto S, Miyano S. Parallel algorithm for learning optimal Bayesian network structure. Journal of Machine Learning Research. 2011;12:2437–2459.
- Tian J, He R. Computing posterior probabilities of structural features in Bayesian networks. Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence; Corvallis, Oregon: AUAI Press; 2009. pp. 538–547.
- Xie X, Geng Z. A recursive method for structural learning of directed acyclic graphs. Journal of Machine Learning Research. 2008;9:459–483.
- Zellner A. On assessing prior distributions and Bayesian regression analysis with g-prior distributions. In: Goel PK, Zellner A, editors. Bayesian Inference and Decision Techniques: Essays in Honour of Bruno de Finetti. North-Holland; Amsterdam: 1986. pp. 233–243.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |