PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of wtpaEurope PMCEurope PMC Funders GroupSubmit a Manuscript
 
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

A Gibbs Sampler for Learning DAGs

Peter Spirtes, Editor

Abstract

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.

Keywords: structure learning, DAGs, Bayesian networks, Gibbs sampling, Markov chain Monte Carlo, variable selection

1. Introduction

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 𝒢 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 MC3 (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).

2. Notation and Model

Let G denote a DAG with vertex set V (G) = {1, … , p}, and directed edge set E(G) [subset or is implied by] 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 AG, with entries specified as AuvG=1(u,v)E(G) and diagonal entries equal to zero. For inference, G is treated as a latent graph in the space 𝒢 of all possible DAGs with p vertices. When a DAG is used to define a probabilistic graphical model (a Bayesian network), each vertex (or node; we use both terms interchangeably) is associated with a component of a p-dimensional random vector X = (X1Xp)T. XZ denotes the random variables (RVs) corresponding to variable indices Z [subset, dbl equals] {1, … , p}. We use bold type for the corresponding data; n samples are collected in the n × p matrix X, with Xu denoting the column corresponding to variable u and XZ the submatrix corresponding to variable index set Z [subset, dbl equals] {1, … , p} (for notational simplicity we assume columns are ordered by increasing variable index).

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 pavG={u:(u,v)E(G)} denote the set of nodes that are parents of node v in G. A tuple of p parent sets pa1=pa1G,,pap=papG fully specifies the edge set E, since upavG(u,v)E(G). Thus, structure learning can be viewed as inference of parent sets pav. The parent set pav takes values in the power set 𝒫(V \ {v}) of nodes excluding node v, subject to acyclicity. We will usually suppress the labels, and simply write pa1G,,papG when specifying parent sets. Since pavG is a subset of the variable indices, we can use XpavG to denote the corresponding data submatrix.

We denote the tuple of parent sets of a set of nodes Z = {z1, … , zs} [subset, dbl equals] V by paZG=(pazG)zZ. This is a tuple of s components (for notational convenience we take these to be ordered by increasing vertex index), each of which is the parent set for the corresponding node in G. We wish to make inferences about paZ, which takes values in 𝒫(V \{z1}) ×…× 𝒫(V \ {zs}), subject to acyclicity. The tuple of parent sets for the complement ZC = V \ Z is denoted paZG=(pazG)zZc. Clearly, a tuple paZ1G,,paZsG of tuples of parent sets specifies the parents of every node in a graph whenever {Z1, … , Zs} forms a partition of V ; the entire edge set can thus be specified by such a tuple. In particular, note that any graph G can be specified as paZG,paZG for any Z [subset or is implied by] V . Some of the notation we use is illustrated in Figure 1.

Figure 1
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 AG as shown, can be specified using parent sets as G = left angle bracketpa1 = [empty], {1}, pa2 = {1}, pa3 ={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(X|G,{θv})=vVp(Xv|XpavG,θv), 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 P(G | X). This is proportional to the product of the marginal likelihood p(X | G) and a graph prior π(G). Our sampler is compatible with essentially any specific model, but inherits computational costs associated with evaluation of the relevant quantities. In all examples we assume conjugate priors for θv, as well as local parameter independence and modularity; this leads to a closed-form marginal likelihood in both multinomial (Heckerman et al., 1995) and Gaussian (Geiger and Heckerman, 1994) cases. Computations are simplified if the graph prior is modular (Friedman and Koller, 2003), meaning it factorises as π(G)=vVπv(pavG); we assume modular priors in all examples. Note that the prior is not specified over the space of orders and so is not subject to the restrictions this entails (Ellis and Wong, 2008; Eaton and Murphy, 2007). Under these assumptions, the posterior factorises across nodes as

P(G|X)=P(pa1=pa1G,,pap=papG|X)vVp(Xv|XpavG)πv(pavG),

where p(Xv|XpavG) is the marginal likelihood for node v given the graph G=pa1G,,papG. The distribution P (G | X) is the target distribution for our sampler.

3. A Gibbs Sampler for Structure Learning

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 MC3) and then describe a naïve Gibbs sampler, which offers no gains over MC3, 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.

3.1. MC3 Sampler

The standard sampler for structure learning of DAGs is MC3 (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,G)=min{1,P(G|X)|neigh(G)|1P(G|X)|neigh(G)|1}.

Variants of MC3 include single-edge direction reversal proposals (Giudici and Castelo, 2003).

3.2. A Naïve (and Inefficient) Gibbs Sampler

Constructing a Gibbs sampler that is analogous to MC3 is straightforward. The posterior distribution on DAGs is a joint distribution over the off-diagonal entries in the adjacency matrix AG, i.e. over the p(p − 1) binary RVs AuvG,uv. At each iteration, MC3 can be thought of as proposing to toggle the value AuvG for some uv, subject to the restriction that the resulting graph must be acyclic.

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 AuvG (for some uv) given the rest of the graph. If (u, v) [set membership] 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) [negated set membership] 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 AuvG is Bernoulli, i.e.

P(AuvG=a|AuvG)={P(G(u,v)|X)P(G(u,v)|X)+P(G(u,v)|X)a=0,P(G(u,v)|X)P(G(u,v)|X)+P(G(u,v)|X)a=1.

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.

3.3. Inefficiency in Gibbs Sampling

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 AuvG, particularly for the collections of RVs corresponding to parent sets. For example, there may be RVs Xu and Xv that in combination score highly as parents of node w, but not individually. Then, AuwG and AvwG will be correlated. In addition, the acyclicity restriction may induce strong dependence. For example, suppose two RVs Xu and Xv are strongly correlated and both edges (u, v) and its reverse (v, u) have high posterior probability. If the edge (u, v) is present, the probability of it being removed is low, but its presence precludes the reversed edge (v, u) from being added. This possibility motivates the edge reversal move used in variants of MC3 (Giudici and Castelo, 2003). Figure 2 shows a three node scenario. Here, both graphs (a) and (b) have high probability. Since reversing the edge (w, x) forms a cycle in (a), moves that consider the parents of only the pair w and x at the same time will not move between graphs (a) and (b) easily. Samplers that alter only a single edge indicator, such as MC3, will also fail. In contrast, sampling the parents of all three nodes jointly would make it easy to move between graphs (a) and (b). Decorrelating transformations could alleviate such problems, but in general finding a suitable transformation is difficult, and for DAGs must of course encapsulate the requirement for acyclicity.

Figure 2
Illustrative scenario in which small local moves limit transitions between two regions of high probability. If both graphs (a) and (b) have high probability, the near-cyclic nature of the graphs makes transitions between (a) and (b) difficult.

3.4. A Gibbs Sampler with Blocking

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 pa1G,,papG of parent sets. It is therefore convenient to describe the sampler using this alternative graph specification, in which G=pa1G,,papG.

We denote the set of q nodes whose parent sets will be sampled together as a block by W = {w1, … , wq} [subset, dbl equals] V. Suppose the current graph is G=paWG,paWG (recall that any graph can be written in this way with respect to any partition of the node set). A move to a new graph G=paWG,paWG is formed by changing the parents of the nodes in W from paWG=(paw1G,,pawqG)topaWG=(paw1G,,pawqG) and setting paWG=paWG (i.e. leaving the parents of nodes not in W unchanged). We sample the tuple paWG of parent sets jointly, conditional on the tuple paWG of parent sets of nodes not in W (that remain unchanged). In terms of the adjacency matrix AG, each block consists of the indicators specifying the parents of the nodes in W, i.e. AvwG for v [set membership] V, w [set membership] W, vw.

To construct a Gibbs sampler using these blocks, we need to find the conditional posterior distribution on the tuple paW of parent sets, given that the remaining tuple paW of parent sets is set to paWG. The conditional distribution depends on whether the graph G′ formed using the proposed parent sets is acyclic. We therefore introduce PaWG to denote the set of permissible tuples, i.e. tuples paWG of parent sets such that G=paWG,paWG is acyclic. For tuples paWGPaWG, the conditional probability is 0. For paWGPaWG, the conditional probability is

P(paW=paWG|paW=paWG,X)=P(paWG,paWG|X)P(paWG|X)=P(G|X)paWGPaWGP(paWG,paWG|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).

Algorithm 1 A Gibbs sampler for learning DAGs, with blocks

Initialise starting point G0=pa1G0,,papG0

for t in 1 to N do

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

    Sample paWGfromP(paW=paWG|paWGt1,X)

    Set GtG=paWG,paWGt1

end for

3.5. Tuning Parameters

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).

3.6. Structure Learning and Variable Selection

It is interesting to note that when q = 1 (and thus W = w for some w [set membership] V) if no choice of parent set induces a cycle then PaWG=𝒫(V\{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 [subset or is implied by] V as parents of node w would create a cycle. Then PaWG=𝒫(V\{w}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.

4. Computational Aspects

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 PaWG 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 PaWG 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.

4.1. A Partition on Permissible Tuples of Parent Sets

We partition the set PaWG of permissible tuples of parent sets as PaWG={PaWG,H1,,PaWG,Hη}. Each component PaWG,Hh 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 Hh [set membership] H, where H is the space of DAGs with q = |W| vertices; η is the cardinality of H. We describe the relationship between each PaWG,Hh and DAG Hh using the following elements, illustrated in Figure 3 and Table 1:

  • The reduced graph G¯=pa1G¯,,papG¯, which is a function of the graph G, and is identical to G except that edges directed into nodes in W are removed.
    pawG¯={wW,pawGwW
  • The (reflexive) transitive closure, which is the directed graph TG on nodes V with edges ET, where (u, v) [set membership] ET 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 z1, … , zs [set membership] V (G), with z1 = u and zs = v, such that (z1, z2), (z2, z3), … , (zs−1, zs) [set membership] E(G)).
  • The descendant nodes dewG¯={vV:TwvG¯=1} and the non-descendant nodes ndwG¯={vV:TwvG¯=0} of w [set membership] W in the reduced graph G¯. Note that wdewG¯andwndwG¯ by definition.
  • The nodes ndG¯=wWndwG¯ that are not descendants in the reduced graph G¯ of any node in W.
  • The nodes dewG¯,H=xpawHdexG¯ that are descendants in the reduced graph G¯ of any node xpawH, for a given node w [set membership] W. Each node xpawH is a parent node of the node w in the graph H=paw1H,,pawqH.
  • The nodes dewG¯,H=xW\pawHdexG¯ that are descenda nts in the reduced graph G¯ of any node xW\pawH, for a given node w [set membership] W. Each node xW\pawH is not a parent node of w in the graph H=paw1H,,pawqH.

Figure 3
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 G¯. Here H is the set of all DAGs on the nodes {3, 4, 7} and thus η = |H| ...
Table 1
The relevant sets and requirements of conditions (A) and (B) for the illustrative graph G shown in Figure 3, with H as shown. Recall ndG¯={1,2}. Condition (B) depends on R7,4G¯,H={4}, but it imposes no restriction when w [set membership] {3, ...

Using this notation, we define PawG,H, for given G [set membership] 𝒢, W [subset, dbl equals] V, w [set membership] W and H=paw1H,,pawqH, as the set of parent sets pawG that satisfy the following conditions.

(A)pawGQwG¯,H=(ndG¯dewG¯,H)\dewG¯,H

(B)pawGRw,xG¯,H for all nodes xpawH, with Rw,xG¯,H=dexG¯\dewG¯,H.

Note that (B) depends on dexG¯ not dexG¯,H. We define PaWG,H as the set of tuples formed by the Cartesian product of the sets PaWG,H of parent sets for w [set membership] W ; in other words paWG=(paw1G,,pawqG)PaWG,H if and only if paw1GPaw1G,H,,pawqGPawqG,H.

Lemma 1 {PaWG,H1,,PaWG,Hη} is a partition of PaWG.

Proof See Appendix B.

Condition (A) ensures all graphs G=paWG,paWG, with parents paWGPaWG,H, are acyclic. It requires that all parents of w [set membership] W are either not descendants in the reduced graph G¯ of any node in W, or a node whose ancestors in the reduced graph 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 PaWG,H for only a single H [set membership] H for any given G [set membership] 𝒢 and thus that {PaWG,H1,,PaWG,Hη} 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 paWG=(paw1G,,pawqG)PaWG,H would not be allowed under any H′ ≠ H, H[set membership] H. Specifically, it ensures that each parent set pawG contains at least one descendant node v [set membership] V of each parent pawH of node w in the graph H, and that v is not a descendant in G¯ of any node x [set membership] W that is not a parent in H of w [set membership] W. To see why this condition is required, consider a graph G′ in which all nodes in W have no parents. Then pawG= for all nodes w [set membership] W and so condition (A) holds for all H [set membership] H. Thus if condition (B) is disregarded, paWGPaWG,H for all H [set membership] H, implying that {PaWG,H1,,PaWG,Hη} is not a partition of PaWG.

4.2. Fast Identification of Partition Components

Parent sets in each set PawG,H can be identified easily using simple set operations, and consequently the partition components PawG,H can be easily identified via Cartesian products of these sets. Let Pawall=𝒫(V\{w}) denote the set of all possible parent sets of node w, subject to maximum in-degree κ, and Pawall(v)=Pawall\𝒫(V\{w,v}) denote the subset all Pawall containing only parent sets that contain node v [set membership] V. Parent sets in PawG,H must (A) not include any nodes not in QwG¯,H, and (B) include at least one node in Rw,xG¯,H for each xpawH, and thus

PawG,H={Paw(B)\Paw(A)ifpawHPawall\Paw(A)ifpawH=,

where Paw(A)=vQwG¯,HPawall(v)andPaw(B)=xpawHrRw,xG¯,HPawall(r).

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 𝒪(pκ+1) time by storing the ‘lookup tables’ Pawall(v) as bit maps, assuming that 𝒪(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 TG with edges ET, the descendants and non-descendants of node u are {v : (u, v) [set membership] ET} and {v : (u, v) [negated set membership] ET} respectively. Thus, the descendants and non-descendants of a node can be identified in 𝒪(1) time when the transitive closure is known. The transitive closure for an arbitrary directed graph with p vertices can be determined in 𝒪(pω) time (Munro, 1971), where ω is the best known exponent for matrix multiplication (Coppersmith and Winograd, 1990, show ω < 2.376). However, in the present setting only incremental changes are made to the current state G of the sampler, so it is not necessary to use an offline algorithm at each iteration. Instead, we can use a fully dynamic transitive closure algorithm (Demetrescu et al., 2010) that provides procedures both for querying the transitive closure, and for incrementally updating it when an edge is added or removed from the graph. We choose to implement the algorithm introduced by King and Sagert (2002), which performs queries in 𝒪(1) time, and updates in 𝒪(p2) worst-case time (see reference for details of required assumptions). A trade-off exists between the performance of these two operations, but this bound for updates is thought to be the best possible whilst retaining 𝒪(1) queries (Demetrescu and Italiano, 2005) and the algorithm is simple to implement.

The algorithm maintains a p × p path count matrix CG whose elements CuvG are the number of distinct paths from node u to node v in the graph G for uv, and CuvG=1 for u = v. The transitive closure can be derived from the path count matrix by noting that TuvG=1 if and only if CuvG>0. Thus query operations are performed in 𝒪(1) by simply checking whether the relevant component of CG is positive. When an edge is added or removed CG is updated as follows. First consider adding an edge (u, v) to a graph G to form a new graph G′. The increase in the number of distinct paths between any two nodes x and y is given by the (x, y) element of C*uGCv*G,whereC*uG denotes the uth column of CG, Cv*G denotes the vth row and [multiply sign in circle] denotes the outer product. Thus, the path count matrix is updated simply as CG=CG+C*uGCv*G. Similarly, when an edge (u, v) is removed from the graph the update is CG=CGC*uGCv*G.

4.3. Two-stage Sampling Method

We draw a new graph G=paWG,paWG starting from a graph G using a two-stage method: we sample first a component PaWG,H of the partition of permissible tuples of parent sets and then a tuple paWG of parent sets from the selected partition component.

In the first stage, PaWG,H is drawn from the conditional distribution, given the tuple of parent sets of nodes not in W.

P(PaWG,H|paWG,X)=P(PaWG,H,paWG|X)P(paWG|X)=paWGPaWG,HwWp(Xw|XpawG)πw(pawG)HpaWGPaWG,HwWp(Xw|XpawG)πw(pawG)=wWpawGPawG,Hp(Xw|XpawG)πw(pawG)HwWpawGPawG,Hp(Xw|XpawG)πw(pawG)

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 [set membership] H, W [subset, dbl equals] V and G [set membership] 𝒢.

paWPaWG,HwWp(Xw|Xpaw)πw(paw)=wWpawPawG,Hp(Xw|Xpaw)πw(paw)

Proof See Appendix C.

In the second stage, we sample new parents paWG from the selected partition component, and form the new graph G=paWG,paWG. The parents of each node w [set membership] W, conditional on H′, are independent, and so can be sampled separately from the following conditional distribution:

P(pawG|PaWG,H,paWG,X)=p(Xw|XpawG)πw(pawG)pawGPawG,Hp(Xw|XpawG)πw(pawG).

This step is straightforward because this distribution is simply the posterior distribution of a constrained variable selection with response w and PawG,H 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 PawG,H, which is 𝒪(pκ+1).

Algorithm 2 A Gibbs sampler for learning DAGs, with general blocks

   Initialise starting point G0=pa1G0,,papG0

   Compute initial path count matrix CG0

   for t in 1 to N do

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

      Let G = Gt−1

      Form G¯ as defined in Section 4.1

      Evaluate CG¯ from CG as described in Section 4.2

      for w [set membership] W do

         Evaluate dewG¯={v:CwvG¯1}

         Evaluate ndwG¯={v:CwvG¯=0}

end for

for H [set membership] H do

   for w [set membership] W do

      Evaluate PawG,H, as described in Section 4.2

      Let KwH=pawGPawG,Hp(Xw|XpawG)πw(pawG)

   end for

   Let KH=wWKwH

end for

Sample PaWG,H according to P(PaWG,H|paWG,X)=KHHKH

for w [set membership] W do

   Sample paWG according to P(pawG|PaWG,H,paWG,X)

end for

Set GtG=paWG,paWGt1

Update CGt

end for

5. Results

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.

5.1. Evaluation Setup

We compare the Gibbs sampler with the MC3 sampler (Madigan and York, 1995) and the REV sampler (Grzegorczyk and Husmeier, 2008), a variant of MC3 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 MC3 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) [proportional, variant] 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 106 iterations of REV (retaining only every 10th iteration to reduce storage requirements). Following Grzegorczyk and Husmeier (2008), 85% of proposals within REV were MC3 proposals (without MC3 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 MC3 moves), but we nevertheless treated the computational costs as the same for the purposes of comparison and drew 106 iterations of the Gibbs sampler (again retaining every 10th iteration). The computational costs of our implementation of MC3 are roughly 1/10 of a Gibbs iteration, and so we performed 107 iterations of MC3 (retaining every 100th).

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.

5.2. Evaluation 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 GMAP, or as a thesholded graph Gτ formed by including all directed edges whose marginal posterior edge probabilities are at least τ (note that Gτ need not be acyclic). When τ = 0.5 we get the median probability model Gmed (for normal linear models Barbieri and Berger, 2004, show that in some settings this is the optimal model for prediction). We also consider thresholding to match sparsity levels seen in frequentist point estimates, i.e. setting τ such that |E(Gτ)|=|E(G^)| for a point estimate G^. By GτPC and GτXie we denote thresholded graphs whose sparsity is matched to the PC and Xie-Geng estimates respectively.

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).

5.3. Small Domain Comparison to Exact Posterior

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 MC3 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 MC3 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 MC3 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 MC3. 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 MC3 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.

Figure 4
For the Zoo data set, maximum and average (across all possible 272 edges) error in posterior edge probabilities versus iteration number (on a log10 scale). For each MCMC algorithm, 10 independent runs, initialised at disparate initial values, are shown. ...

5.4. The ALARM Network

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 MC3 −11,176.9 ± 273.7. The highest scoring graph discovered by any of the 10 runs of the MC3 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.

Figure 5
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 log10 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 MC3 and REV respectively (when n = 5000).

Figure 6
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 MC3 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 MC3 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 MC3 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).

Figure 7
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, ...
Figure 8
ALARM data, area under the receiver operating characteristic curve (AUROC) against iteration number (log10 scale) for each of the 10 independent runs for each MCMC algorithm. Each panel shows results for a particular sample size n = 100, … , 5000. ...

5.5. Networks of Varying Sparseness

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 AG, by drawing entries as AuvGBernoulli(ρ), where ρ [set membership] (0, 1) is a parameter controlling sparsity. Entries were drawn independently for each u < v, with AuvG=0 otherwise. Larger ρ gives a denser graph and the expected number of neighbours (parents or immediate children) of a node is ρ(p−1). We set p = 25, n = 1000 and considered 5 values of the sparseness parameter ρ (corresponding to expected neighbourhood sizes 2, 3, 4, 5 and 6). For each ρ we drew 10 DAGs, simulating a data set from each DAG by ancestral sampling using a Normal linear model (see Kalisch and Bühlmann, 2007, for full details).

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 Gmed 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.

Figure 9
Synthetic data, accuracy of estimation for networks with different levels of sparsity. The panels correspond to simulations in which the expected number of neighbours for each node in the data-generating graph is respectively 2, 3, 4, 5, and 6. Accuracy ...

5.6. Survey Data

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 Gmed 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 MC3 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 MC3. The Gibbs sampler also had the highest proportion of edges with PSRF < 1.1 (Figure A14c).

Figure 10
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 103 samples were needed). REV took longer to reach a (usually lower) plateau. MC3 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 MC3. For example, using Gτ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 (MC3) and 30.8 (PC-algorithm). Results for GMAP and GτXie are shown in Figure A19a.

5.7. Large-sample, Single-cell Molecular Data

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 Gmed 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 MC3. 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 MC3 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 (MC3). 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τPC point estimates) were 128.8 (Gibbs), 140.1 (REV) and 225.2 (MC3). We note that the SHDs are particularly high here because the PC-algorithm estimates a network with a high density (recall that under GτPC all methods choose the same number of edges as PC).

Figure 11
Convergence diagnostic plots for each MCMC sampler for the single-cell molecular data. The posterior edge probabilities given by two independent MCMC runs are plotted against each other, binned into hexagonal areas to avoid overplotting. When the two ...

5.8. A Pathological 4-node Example

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 105 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 106 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 106 iterations), but in real-world examples it is difficult to rule out the possibility that analogous issues may arise.

Figure 12
For the pathological 4-node example, maximum and average (across all possible 12 edges) error in posterior edge probabilities by iteration number (on a log10 scale). For each MCMC algorithm, 10 independent runs, initialised at disparate initial values, ...

6. Discussion

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(paW|paWG,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(paW|paWG,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 MC3. 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 MC3 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 MC3 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(paW|paWG,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(paW|paWG,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, MC3 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.

Supplementary Material

Appendix

Acknowledgments

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).

Contributor Information

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

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

References

  • 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(n2) 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.