Science, Faculty of
Computer Science, Department of
DSpace
UBCV
Eaton, Daniel James
Master of Science - MSc
University of British Columbia
In this work, we address both the computational and modeling aspects of Bayesian network structure learning. Several recent algorithms can handle large networks by operating on the space of variable orderings, but for technical reasons they cannot compute many interesting structural features and require the use of a restrictive prior. We introduce a novel MCMC method that utilizes the deterministic output of the exact structure learning algorithm of Koivisto and Sood to construct a fast-mixing proposal on the space of DAGs. We show that in addition to fixing the order-space algorithms' shortcomings, our method outperforms other existing samplers on real datasets by delivering more accurate structure and higher predictive likelihoods in less compute time. Next, we discuss current models of intervention and propose a novel approach named the uncertain intervention model, whereby the targets of an intervention can be learned in parallel to the graph's causal structure. We validate our model experimentally using synthetic data generated from known ground truth. We then apply our model to two biological datasets that have been previously analyzed using Bayesian networks. On the T-cell dataset of Sachs et al. we show that the uncertain intervention model is able to better model the density of the data compared to previous techniques, while on the ALL dataset of Yeoh et al. we demonstrate that our method can be used to directly estimate the genetic effects of the disease.
Bayesian network structure learning for the uncertain experimentalist W i t h applications to network biology by Daniel James Eaton B.Sc, The University of British Columbia, 2005 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF Master of Science in The Faculty of Graduate Studies (Computer Science) The University Of British Columbia June, 2007 © Daniel James Eaton 2007 Abstract In this work, we address both the computational and modeling aspects of Bayesian network structure learning. Several recent algorithms can handle large networks by operating on the space of variable orderings, but for tech-nical reasons they cannot compute many interesting structural features and require the use of a restrictive prior. We introduce a novel M C M C method that utilizes the deterministic output of the exact structure learning algorithm of Koivisto and Sood to construct a fast-mixing proposal on the space of DAGs. We show that in addition to fixing the order-space algorithms' shortcomings, our method outperforms other existing samplers on real datasets by delivering more accurate structure and higher predictive likelihoods in less compute time. Next, we discuss current models of intervention and propose a novel approach named the uncertain intervention model, whereby the targets of an interven-tion can be learned in parallel to the graph's causal structure. We validate our model experimentally using synthetic data generated from known ground truth. We then apply our model to two biological datasets that have been previously analyzed using Bayesian networks. On the T-cell dataset of Sachs et al. we show that the uncertain intervention model is able to better model the density of the data compared to previous techniques, while on the A L L dataset of Yeoh et al. we demonstrate that our method can be used to directly estimate the genetic effects of the disease. ii Table of Contents Abstract 1 1 Table of Contents n i List of Tables v List of Figures v l Acknowledgements v n l Dedication l x I Thesis 1 1 Introduction 2 1.1 Motivation 2 1.2 Outline of Thesis 4 2 Structure learning by dynamic programming and M C M C . . 6 2.1 Introduction 6 2.2 Previous work 8 2.3 Modular priors 10 2.4 Our method 13 2.4.1 Likelihood models 15 2.5 Experimental results 18 2.5.1 Speed of convergence to the exact posterior marginals . . 18 iii Table of Contents 2.5.2 Structural discovery 19 2.5.3 Accuracy of predictive density 20 2.5.4 Convergence diagnostics 26 2.6 Summary and future work 26 3 Structure learning with uncertain interventions 27 3.1 Introduction 27 3.2 Models of intervention 28 3.2.1 No interventions 29 3.2.2 Perfect interventions . . . 30 3.2.3 Imperfect interventions 31 3.2.4 Unreliable interventions 32 3.2.5 Soft interventions 33 3.2.6 Uncertain interventions 34 3.2.7 The power of interventions 36 3.3 Algorithms for structure learning 38 3.3.1 Exact algorithms for ;;(HY,-) and H * M A P 38 3.3.2 Local search for UMAP 39 3.3.3 Iterative algorithm for HMAP 39 3.4 Experimental results 40 3.4.1 Synthetic data 40 3.4.2 T-cell data 46 3.4.3 A L L data 53 3.5 Summary and future work 59 Bibliography 63 II A p p e n d i c e s 69 A Software 70 iv List of Tables 1.1 Number of DAGs on d nodes 3 A . l Major features of BNSL software 70 v List of Figures 2.1 Deviation of various modular priors from uniform 12 2.2 Cancer network 15 2.3 Convergence to edge marginals on Cancer network 16 2.4 Convergence to edge marginals on CHD dataset 16 2.5 Edge recovery performance of M C M C methods on Child network 17 2.6 Path recovery performance of M C M C methods on Child network 17 2.7 Child network 20 2.8 Test set log likelihood vs training time on Synthetic-15 network . 21 2.9 Test set log likelihood vs training time on Adult dataset 21 2.10 Test set log likelihood vs training time on T-cell dataset 22 2.11 Samplers' training set log likelihood trace plots on Adult dataset 25 3.1 Intervention model: None 29 3.2 Intervention model: Perfect 30 3.3 Intervention model: Imperfect 31 3.4 Intervention model: Imperfect with unreliable extension 32 3.5 Intervention model: Soft 33 3.6 Example of a "Fat hand" intervention 36 3.7 Markov equivalence classes on the Cancer network 37 3.8 Perfect vs uncertain interventions on the Cancer network . . . . 42 3.9 Cars network ground truth 43 3.10 Structural recovery performance of perfect vs uncertain interven-tion models on Cars network 45 vi List of Figures 3.11 Structural recovery performance of exact B M A vs point estima-tion on Cars network 47 3.12 T-cell dataset 48 3.13 T-cell dataset results 50 3.14 Predictive likelihood scores of intervention models on T-cell dataset 51 3.15 A L L dataset 54 3.16 Highest scoring DAG found on A L L dataset 55 3.17 Predicted expression profile for a patient without A L L 56 3.18 Sampled expression profile for a patient with A L L subtype E2A-PBX1, M L L or both 57 3.19 A L L subtypes E2A-PBX1 and B C R - A B L results 60 3.20 A L L subtypes Hyperdip > 50, TEL-AML1, M L L and T - A L L results 61 3.21 Other A L L subtypes' results 62 Acknowledgements I owe at least two sigmas of thank you density to my supervisor Dr. Kevin Murphy. Over the past two years he has been an incredible source of inspiration and motivation, keeping the Bayes ball rolling even in the face of occasionally disappointing results or formidable bugs. I found the support of an NSERC scholarship to be most beneficial, and therefore thank Canada for supporting my work. I also appreciate the instruction of several professors in the Computer Sci-ence department, and the help I received from peers while preparing papers for conference submission, especially Mark Schmidt, Hoyt Koepke and Hendrik Kueck. I would like to acknowledge my family's contributions: Joan valiantly proof-read my thesis, while Sheena, Doug and Sarah provided me with motivation to graduate on time, some in more gentle ways then others. I also wish to thank Michael Forbes for generously sharing his UBC thesis I^TpX template; it made formatting a breeze. Daniel Eaton University of British Columbia June 2007 vm Dedica t ion This thesis is dedicated to my parents, Doug and Joan, in recognition of their continuous support through each episode of my academic career. ix Part I Thesis Chapter 1 Introduction 1.1 Motivation Bayesian networks are a powerful tool for modeling stochastic systems, as they enable high-dimensional joint probability distributions to be represented effi-ciently The model is especially popular due to the availability of off-the-shelf algorithms and software to learn them from data, and then to extract mean-ingful quantities in support of scientific hypotheses. More importantly, unlike undirected models, the inter-variable relationships encoded by a Bayesian net-work can be interpreted causally, meaning that the notion of "smoking causes cancer" can be characterized directly, rather than simply, "smoking is linked to cancer". This distinction may appear subtle, but the notion of causality is at the heart of modern science [45]. Therefore, it comes as no surprise that Bayesian networks are rapidly spreading beyond the domains for which they were originally created. The systems biology community has recently begun exploiting the model to analyze high-throughput biological data, such as gene expression microarrays [12, 20, 25]. In many classical applications, the relationships between variables were assumed known a priori; in other words, the Bayesian network structure was given, perhaps by a human expert. The goal then became to learn the associated model parameters. In contrast, many system biologists aim to learn the structure of the network from data, for example to infer gene regulatory in-teractions, discover new drugs, or to guide costly future experiments. Evidence suggests that cancers and many other diseases manifest themselves by disrupting specific cellular processes [1]. Gene expression microarrays can measure tran-2 Chapter 1. Introduction scription levels for thousands of genes in a cell, offering insight into the complex cell dynamics that give rise to diseases. Causal Bayesian network models supply a natural means of modeling systems biology data, potentially yielding answers about cellular processes, or identify-ing the effects of external agents, such as diseases or drugs. In order to uniquely infer causal relationships, it is necessary to perform experiments where variables are perturbed or manipulated to particular known values. In biological appli-cations, an intervention might be accomplished by injecting chemicals into a cell or performing a specific gene knockout. Alternatively, nature provides un-wanted perturbations in the form of disease, which can be construed as virtual experiments. Unfortunately, even if both observational and experimental data is readily available, Bayesian network structure learning remains a challenging problem, largely due to the super-exponential growth of the space of valid networks given the number of variables. Valid structures are represented by directed acyclic graphs (DAGs), and there are O(d!20) DAGs on d nodes [47], Table 1.1 shows the growth for small d. Exhaustive search or naive marginalization over such a space is clearly intractable in general. Therefore, without clever algorithms it is necessary to resort to local search through a huge hypothesis space. d 2 3 4 5 6 7 8 9 #DAG(d) 3 25 543 29281 3781503 l.le9 7.8ell 1.2el5 Table 1.1: Number of DAGs on d nodes The scarcity of data inherent to most biological applications greatly exacer-bates this computational problem. Often, the number of variables dwarfs the sample size, meaning that many DAGs are likely to fit the data well. In this setting, it would be dangerous to commit to a particular structure and make any interpretation on the causal relationships between variables, since there conceivably could be another structure that also fits the data well, and leads to contradictory conclusions. Full Bayesian model averaging would be the greatly 3 Chapter 1. Introduction preferable strategy; however, the naive approach becomes intractable for more than a mere 6 variables. Lastly, many of the scientific questions that the systems biology community poses are unanswerable using existing models of intervention. Typically when a system is intervened on, the experimenter assumes the targets of their pertur-bation are known, and exploits this knowledge in the prescribed way [8] during structure learning. The perturbation may instead have a "fat hand" and cause unexpected side-effects that cannot be explained by the intrinsic relationships between variables. Perhaps the goal is to learn the intervention's targets, for example to measure the side-effects of a new treatment. Unfortunately, it is not possible to directly answer these, nor other important questions with existing models. 1.2 Outline of Thesis This thesis presents novel algorithms and models for Bayesian network structure learning aimed at overcoming the challenges introduced in Section 1.1. In Chap-ter 2, we discuss existing computational methods, including Markov chain monte carlo (MCMC) and exact Bayesian model averaging (BMA), and briefly discuss their shortcomings. Next, we introduce an algorithm that combines M C M C with exact B M A to solve them, and show that it simultaneously outperforms the other sampling methods based on time and accuracy. Chapter 3 presents a new model of experimental data, allowing for uncertainty in the targets of interventions. The new model is verified on synthetic data, and then tested on two gene expression datasets, producing results that agree with other reported analyses, but also providing novel benefits that the other methods cannot. The chapter also shows how to adapt a recent exact structure learning algorithm to handle experimental data. Chapters 2 (based on [13]) and 3 (based on [14]) are structured in a self-contained format, each containing background, results and conclusions relevant to their respective topics. However, it must be noted that each topic is heavily 4 Chapter 1. Introduction intertwined: algorithms borrowing models to test on data, and models using algorithms to do computation. 5 Chapter 2 Structure learning by dynamic programming and 2.1 Introduction Directed graphical models are useful for a variety of tasks, ranging from density estimation to scientific discovery. One of the key challenges is to learn the structure of these models from data. Often (e.g., in molecular biology) the sample size is quite small relative to the size of the hypothesis space. In such cases, the posterior over graph structures given data, p(G\D), gives support to many possible models, and using a point estimate (such as MAP) could lead to unwarranted conclusions about the structure, as well as poor predictions about future data. It is therefore preferable to use Bayesian model averaging. If we are interested in the probability of some structural feature / (e.g., f{G) = 1 if there is an edge from node i to j and /(G) = 0 otherwise), we can compute the posterior mean estimate E(f\D) = J2G f{G)p(G\D). Similarly, to predict future data, we can compute the posterior predictive distribution p(x\D) = exact Bayesian model averaging is intractable in general. However, if the model space is restricted special cases arise where full averaging is possible; for example, over all trees [40] or all graphs consistent with a given node ordering [10]. If M C M C Y:GP(X\G)P(G\D). DAGs (directed acyclic graphs) on d nodes [47], 6 Chapter 2. Structure learning by dynamic programming and MCMC the ordering is unknown, we can use M C M C techniques to sample orders, and sample DAGs given each such order [19, 21, 29]. However, Koivisto and Sood [31, 32] showed that one can use dynamic programming (DP) to marginalize over orders analytically. This technique enables one to compute all marginal posterior edge probabilities, p(Gij = 11-0), exactly in 0(d2d) time. Although exponential in d, this technique is quite practical for d < 20, and is much faster than comparable M C M C algorithms on similar sized problems1. Unfortunately, the DP method has three fundamental limitations, even for small domains. The first problem is that it can only be used with certain kinds of graph priors which satisfy a "modularity" condition, which will be described in Section 2.3. Although this seems like a minor technical problem, it can result in significant bias. This can lead to unwarranted conclusions about structure as well as poor predictive performance, even in the large sample setting. The second problem is that it can only compute posteriors over modular features; thus it cannot be used to compute the probability of features like "is there a path between nodes i and j via fc", or "is i an ancestor of j". Such long-distance features are often of more interest than direct edges. The third problem is that it is expensive to compute predictive densities, p(x\D). Since the DP method integrates out the graph structures, it has to keep all the training data D around, and predict using p(x\D) = p(x, D)/p(D). Both terms can be computed exactly using DP, but this requires re-running DP for each new test case x. In addition, since the DP algorithm assumes complete data, if x is incompletely observed (e.g., we want to "fill in" some of it), we must run the DP algorithm potentially an exponential number of times. For the same reason, we cannot sample from p(x\D) using the DP method. We propose to fix all three of these shortcomings by combining DP with the Metropolis-Hastings (MH) algorithm. The basic idea is simply to use the 1 Our Mat lab /C implementation takes 1 second for d = 10 nodes and 6 minutes for d = 20 nodes on a standard laptop. The cost is dominated by the marginal likelihood computation, which all algorithms must perform. Our code is freely available; please see Appendix A for instructions to obtain it. 7 Chapter 2. Structure learning by dynamic programming and MCMC DP algorithm as an informative (data driven) proposal distribution for mov-ing through DAG space, thereby getting the best of both worlds: a fast de-terministic approximation, plus unbiased samples from the correct posterior, Gs ~ p(G\D). These samples can then be used to compute the posterior mean of arbitrary features, E[f\D] « ^ Z)f=i f(Gs)> o r t n e posterior predictive distri-bution, p(x\D) « ^ J2s=i p(a^|Cs). Results presented in Section 2.5 show that this hybrid method produces more accurate estimates than other approaches, given a comparable amount of compute time. The idea of using deterministic algorithms as a proposal has been explored before e.g. [11], but not, as far as we know, in the context of graphical model structure learning. Further, in contrast to [11], our proposal is based on an exact algorithm rather than an approximate algorithm. 2.2 Previous work The most common approach to estimating the posterior p(G\D), or marginal features thereof, is to use the Metropolis-Hastings (MH) algorithm, using a proposal that randomly adds, deletes or reverses an edge; this has been called MC3 for Markov Chain Monte Carlo Model Composition [36]. (See also [24] for some improvements, and [35] for a related approach called Occam's window.) Unfortunately, this proposal is very local, and the resulting chains do not mix well in more than about 10 dimensions. An alternative is to use Gibbs sampling on the adjacency matrix [37], In our experience, this gets "stuck" even more easily, although this can be ameliorated somewhat by using multiple restarts, as the experimental results will demonstrate. A different approach, first proposed in [21], is to sample in the space of node orderings using MH. It is based on the fact that conditioned on a node ordering -<, the probability of the data and a feature factorizes d P(XJ\I,^) = H J2 p ( G l | ^ ) p ( X 1 | X G , , / ) / I ( G I ) (2.1) 8 Chapter 2. Structure learning by dynamic programming and MCMC where p(x,f) means p(x,f(G) = 1) and — {j : j -< i] are the set of nodes that preceed i. This fact was first observed by [3], and was exploited by [21] using a M H proposal that randomly swaps the ordering of nodes. For example, (1,2,3,4,5,6) -» (1,5,3,4,2,6) where we swapped 2 and 5. This is a smaller space ("only" 0(d\)), and is "smoother", allowing chains to mix more easily. [21] provides experimental evidence that this approach gives much better results than M H in the space of DAGs with the standard add/ delete/ reverse proposal. Unfortunately, in order to use this method, one is forced to use a modular prior, which has var-ious undesirable consequences that we discuss in Section 2.3. Ellis and Wong [19] realized this, and suggested using an importance sampling correction. How-ever, computing the exact correction term is #P-hard, and our empirical results suggest that their approximate correction yields inferior structure learning and predictive density compared to our method. An alternative to sampling orders is to analytically integrate them out using dynamic programming (DP) [31, 32]. The algorithm is complex, but the key idea can be stated simply: when considering different variable orderings — say (3,2,1) and (2,3,1) — the contribution to the marginal likelihood for some nodes can be re-used. For example, p(Xi\X2, X:i) is the same as p(X-l\X3,X2), since the order of the parents does not matter. By appropriately caching terms, one can devise an 0(d2d) algorithm to exactly compute the marginal likelihood and marginal posterior features. The inputs to this algorithm are a modular prior (see Section 2.3) and the local conditional marginal likelihoods, which must be computed for every node i and every possible parent set Gj (up. to size k). There are y_p_0 (^ ) = 0(dk) such terms, therefore the time complexity of the DP algorithm including marginal likelihood computation is 0(d2d + dk+lC(N)). Here, C(N) is the amount of time needed to compute each marginal likelihood term as a function of the sample size N. In practise, the polynomial term usually strongly dominates the exponential term. To compute the posterior predictive density, p(x\D), the standard approach 9 Chapter 2. Structure learning by dynamic programming and MCMC is to use a plug-in estimate p(x\D) ~ p(x\G(D)). Here G may be an approximate M A P estimate computed using local search [27], or the MAP-optimal DAG which can be found by the recent algorithm of [50] (which takes o(d22d~2) time.) Alternatively, G could be a tree; this is a popular choice for density estimation since one can compute the optimal tree structure in 0(d2 logd) time [7, 41]. It can be proven that averaging over the uncertainty in G will, on average, produce higher test-set predictive likelihoods [34]. The DP algorithm can com-pute the marginal likelihood of the data, p(D) (marginalizing over all DAGs), and hence can compute p{x\D) = p(x,D)/p(D) by calling the algorithm twice. (We only need the "forwards pass" of [32], using the feature / = 1; we do not need the backwards pass of [31].) However, this is very expensive, since we need to compute the local marginal likelihoods for every possible family on the expanded data set for every test case x. In Section 2.5 we will show that by averaging over a sample of graphs our method gives comparable predictive performance at a much lower cost. 2.3 Modular priors Some of the best current methods for Bayesian structure learning operate in the space of node orders rather than the space of DAGs, either using M C M C [19, 21, 29] or dynamic programming [31, 32]. Rather than being able to define an arbitrary prior on graph structures p(G), methods that work with orderings define a joint prior over graphs G and orders -<! as follows: 1 d p(*,G) = zY[qi{U?)Pi(Gi) x /(consistent^, G)) t=i where is the set of predecessors (possible parents) for node i in -<, and Gi is the set of actual parents for node i. We say that a graph structure G = ( G i , . . . , G,j) is consistent with an order (Ui,..., £/,/) if Gj C Ui for all i. (In addition we require that G be acyclic, so that -< exists.) Note that Ui and Gi are not independent. Thus the <?* and pi terms can be thought of as factors 10 Chapter 2. Structure learning by dynamic programming and MCMC or constraints, which define the joint prior p(-<,G). This is called a modular prior, since it decomposes into a product of local terms. It is important for computational reasons that Pi(Gj) only give the prior weight to sets of parents, and not to their relative order, which is determined by Qi(Ui). This feature is what enables the order space algorithms to re-use scores for all orderings of a parent set, and turn a sum over the super-exponential structure space into a sum over the factorial order space, and then further reduce it to an exponential complexity. From the joint prior, we can infer the marginal prior over graphs, p(G) = y__,p(-<1G). Unfortunately, this prior favors graphs that are consistent with more orderings. For example, the fully disconnected graph is the most probable under a modular prior, and trees are more probable than chains, even if they are Markov equivalent (e.g., 1<—2—>3 is more probable than 1—>2—>3). This can cause problems for structural discovery. To see this, suppose the sample size is very large, so the posterior concentrates its mass on a single Markov equivalence class. Unfortunately, the effects of the prior are not "washed out", since all graphs with the equivalence class have the same likelihood. Thus we may end up predicting that certain edges are present due to artifacts of our prior, which was merely chosen for technical convenience. In the absence of prior knowledge, one may want to use a uniform prior over DAGs 2 . However, this cannot be encoded as a modular prior. To see this, let us use a uniform prior over orderings, Qi{Ui) = 1, so p(^) = l/(d\). This is reasonable since typically we do not have prior knowledge on the order. For the parent factors, let us use pi(Gi) = 1; we call this a "modular flat" prior. However, this combination is not uniform over DAGs after we sum over orderings: see Figure 2.1. A more popular alternative (used in [19, 21, 31, 32]) is to take /9j(Gj) oc (j^.j) we call this the "Koivisto" prior. This prior says that different cardinalities of parents are considered to be equally likely a priori. 2 One could argue that we should use a uniform prior over P D A G s , but we will often be concerned with learning causal models from interventional data, in which case we have to use DAGs. 11 Chapter 2. Structure learning by dynamic programming and MCMC Koivisto: KL from uniform = 2.82 Figure 2.1: Deviation of various modular priors from uniform. Priors on all 29,281 DAGs on 5 nodes. Koivisto prior means using Pi(Gi) oc • Ellis prior means the same pit but dividing by the number of consistent orderings for each graph (computed exactly). Modular flat means using p,(Gi) oc 1. This is the closest to uniform in terms of K L distance, but will still introduce artifacts. If we use Pi(Gi) oc 1 and divide by the number of consistent orders, we will get a uniform distribution, but computing the number of consistent orderings is #P-hard. 12 Chapter 2. Structure learning by dynamic programming and MCMC However, the resulting p(G) is even further from uniform: see Figure 2.1. Ellis and Wong [19] recognized this problem, and tried to fix it as follows. Let p*{G) = ^ Yli Pi{G) o e the desired prior, and let p{G) be the actual modular prior implied by using p* and g, = 1. We can correct for the bias by using an importance sampling weight given by w(G) P*{G) jII.PtiGi) P(G) i Ui P*(G i)/(consistent(^, G)) If we set p* = 1 (the modular flat prior), then this becomes w(G) = •=——, : ', E ^ / ( c O n S l S t e n t ( - < , G ) Thus this weighting term compensates for overcounting certain graphs, and in-duces a globally uniform prior, p(G) oc 1. However, computing the denom-inator (the number of orders consistent with a graph) is #P-complete [2]. Ellis and Wong approximated this sum using the sampled orders, w(G) m i^r————. However, these samples ~<s are drawn from the poste-rior p(-< \D), rather than the space of all orders, so this is not an unbiased esti-mate. Also, they used Pi(Gi) cc (j'j.j) \ rather than p{ = 1, which still results in a highly non-uniform prior, even after exact reweighting (see Figure 2.1). In contrast, our method can cheaply generate samples from an arbitrary prior. 2.4 Our method As mentioned above, our method is to use the Metropolis-Hastings algorithm with a proposal distribution that is a mixture of the standard local proposal, that adds, deletes or reverses an edge at random, and a more global proposal that uses the output of the DP algorithm: _ J qiocai(G'\G) w.p. B [ q<jlobal{G ) W .p . l - B . The local proposal chooses uniformly at random from all legal single edge additions, deletions and reversals. Let nbd(G) denote the set of acyclic neighbors 13 Chapter 2. Structure learning by dynamic programming and MCMC generated in this way. We have the proposal distribution The global proposal includes an edge between i and j with probability ptj + Pji < 1, where pij = p{Gij\D) are the exact marginal posteriors computed using DP (using a modular prior). If this edge is included, it is oriented as i—>j w.p. qtj = pij/(pij +Pji), otherwise it is oriented as z<— j . After sampling each edge pair, we check if the resulting graph is acyclic. (The acyclicity check can be done in amortized constant time using the ancestor matrix trick [24].) This leads to If we set (3 — 1, we get the standard local proposal. If we set ft = 0, we get a purely global proposal. Note that qgi0bai{G') is independent of G, so this is an independence sampler. We tried various other settings of ft (including adapting it according to a fixed schedule), which results in performance somewhere in between purely local and purely global. For /3 > 0 the chain is aperiodic and irreducible, since the local proposal has both properties [36, 46]. However, if /? = 0, the chain is not necessarily aperiodic and irreducibile, since the global proposal may set py = pji = 0. This problem is easily solved by truncating edge marginals which are too close to 0 or 1, and making the appropriate changes to qgioi>ai{G'). Specifically, any p^ < C is set to C, while p^ > 1 — C are set to 1 — C. We used C = le — 4 in our experiments. We then accept moves proposed from qgi0bai{G') with probability 14 Chapter 2. Structure learning by dynamic programming and MCMC Figure 2.2: Cancer network. First reported in [22], 2.4.1 Likelihood models For simplicity, we assume all the conditional probability distributions (CPDs) are multinomials (tables), p(Xi = k\Xci = 3,8) — Qijk, though our software (see Appendix A) can handle the linear-Gaussian case as well. We make the usual assumptions of parameter independence and modularity [27], and we use uniform conjugate Dirichlet priors #y ~ Dir(cti,..., a;), where we set a* = l/(9iri)> where qi is the number of states for node Xj and rt is the number of states for the parents Xd- The resulting marginal likelihood, can be computed in closed form, and is called the BDeu (Bayesian Dirichlet likelihood equivalent uniform) score [27]. We use AD trees [43] to compute these terms efficiently. Note that our technique can easily be extended to other CPDs (e.g., decision trees [5]), provided p(Xi\XGi) can be computed or approximated (e.g., using BIC). p ( D \ G ) = \ \ p i X i \ X c , ) 15 Chapter 2. Structure learning by dynamic programming and MCMC Cancer s « e •a * i iS 0 0 30 40 60 80 100 120 HO 0.8 Q <U 0.4 0 20 40 60 80 100 120 140 T ime (seconds) Figure 2.3: Convergence to edge marginals on Cancer network. SAD error vs running time on the 5 node Cancer network (shown in Figure 2.2). The Gibbs sampler performs poorly, therefore we replot the graph with it removed (bottom figure). Note that 140 seconds corresponds to about 130,000 samples from the hybrid sampler. The error bars (representing one standard deviation across 25 chains starting in different conditions) are initially large, because the chains have not burned in. This figure is best viewed in colour. Local Global Hybrid Coronary 0 50 100 150 200 T i m e (seconds) Figure 2.4: Convergence to edge marginals on C H D dataset. Similar to Figure 2.3, but on the 6 node CHD dataset. L6 Chapter 2. Structure learning by dynamic programming and MCMC Edge Features ? 0.85 Local Order Global Hybrid Figure 2.5: Edge recovery performance of M C M C methods on Child network. Area under the ROC curve (averaged over 10 M C M C runs) for detecting edge presence for different methods on the d = 20 node Child network with n = 10/c samples using 200 seconds of compute time. The AUC of the exact DP algorithm is indistinguishable from the global method and hence is not shown. Path Features 0.9 3 < 0.8 0.7 Figure 2.6: Path recovery performance of M C M C methods on Child network. Area under the ROC curve (averaged over 10 M C M C runs) for detecting path presence for different methods on the d = 20 node Child network with n = 10/c samples using 200 seconds of compute time. 17 Chapter 2. Structure learning by dynamic programming and MCMC 2.5 Experimental results 2.5.1 Speed of convergence to the exact posterior marginals In this section we compare the accuracy of different algorithms in estimating p(Gij = 1\D) as a function of their running time, where we use a uniform graph prior p(G) oc 1. (Obviously we could use any other prior or feature of interest in order to assess convergence speed, but this seemed like a natural choice, and enables us to compare to the raw output of DP.) Specifically, we compute the sum of absolute differences (SAD), St = J2ij = 1|-D) - qt{Gi:i = l\D)\, versus running time t, where p(G;j = 1|D) are the exact posterior edge marginals (computed using brute force enumeration over all DAGs) and qt(Gij\D) is the approximation based on samples up to time t. We compare 5 M C M C methods: Gibbs sampling on elements of the adjacency matrix, purely local moves through DAG space (8 = 1), purely, global moves through DAG space using the DP proposal (8 = 0, which is an independence sampler), a mixture of local and global (probability of local move is 8 = 0.1), and an M C M C order sampler [21] with Ellis' importance weighting term.3 (In the figures, these are called as follows: 8 = 1 is "Local", 8 = 0 is "Global", 8 = 0.1 is "Hybrid".) In our implementation of the order sampler, we took care to implement the various caching schemes described in [21], to ensure a fair comparison. However, we did not use the sparse candidate algorithm or any other form of pruning. For our first experiment, we sampled data from the 5 node "Cancer network" of [22] (shown in Figure 2.2) and then ran the different methods. In Figure 2.3, we see that the DP+MCMC samplers outperform the other samplers. We also ran each method on the well-studied coronary heart disease (CHD) dataset [18]. This consists of about 200 cases of 6 binary variables, encoding such things as "is your blood pressure high?", "do you smoke?", etc. In Figure 2.4, we see 3 Without tlie reweighting term, tlie M C M C order sampler [21] would eventually converge to the same results (as measured by SAD) as the DP method [31, 32], 18 Chapter 2. Structure learning by dynamic programming and MCMC again that our DP+MCMC method is the fastest and the most accurate. 2.5.2 Structural discovery In order to assess the scalability of our algorithm, we next looked at data gen-erated from the 20 node "Child" network used in [55] (shown in Figure 2.7). We sampled n = 10, 000 records using random multinomial CPDs sampled from a Dirichlet, with hyper-parameters chosen by the method of [6], which ensures strong dependencies between the nodes4. Stronger dependencies increase the likelihood that the distribution will be numerically faithful to the conditional independency assumptions encoded by the structure. Next, we compute the posterior over two kinds of features: edge features, fij = 1 if there is an edge between i and j (in either orientation), and path features, fij = 1 if there is a directed path from i to j. (Note that the latter cannot be computed by DP; to compute it using the order sampler of [21] requires sampling DAGs given an or-der.) We can no longer compare the estimated posteriors to the exact posteriors (since d = 20), but we can compare them to the ground truth values from the generating network. Following [28, 31], we threshold these posterior features at different levels, to trade off sensitivity and specificity We summarize the resulting ROC curves in a single number, namely area under the curve (AUC). The results for edge features are shown in Figure 2.5. We see that the D P + M C M C methods do very well at recovering the true undirected skeleton of the graph, obtaining an AUC of 1.0 (same as the exact DP method). We see that our D P + M C M C samplers are significantly better (at the 5% level) than the DAG sampler and the order sampler. The order sampler does not do as well as the others, for the same amount of run time, since each sample is more expensive to generate. The results for path features are shown in Figure 2.6. Again we see that 4 For example, consider a node I with 3 states and 4 parent states. The method of [6] prescribes that we pick a "basis vector" (1,1/2,1/3), and then, for the j ' t h parent state, we sample ~ Dir(sctij.), where a u . cc (1,1/2,1/3), a i 2 - oc (1 /2,1 /3,1), a i 3 . oc (1/3,1, 1/2), an. oc (1,1 /2,1 /3), and s = 10 is an effective sample size. 19 Chapter 2. Structure learning by dynamic programming and MCMC Figure 2.7: Chi ld network. First reported in [9]. the D P + M C M C method (using either B = 0 or B — 0.1) yields statistically significant improvement (at the 5% level) in the AUC score over other M C M C methods on this much harder problem. 2.5.3 A c c u r a c y of predict ive density In this section, we compare the different methods in terms of the log loss on a test set: .. rn i=E\ogp(x\D) * - V l o g p f o l D ) m '—' where m is the size of the test set and D is the training set. This is the ultimate objective test of any density estimation technique, and can be applied to any dataset, even if the "ground truth" structure is not known. The hypothesis that we wish to test is that methods which estimate the posterior p(G\D) more accurately will also perform better in terms of prediction. We test this hypoth-esis on three datasets: synthetic data from a 15-node network, the "Adult" US census dataset from the UC Irvine repository and a biological dataset related to the human T-cell signalling pathway [48]. 20 Chapter 2. Structure learning by dynamic programming and MCMC Synthetic-15 -0.25 •0.3 — — - — • i.—I—i h — t - - I • Glbbs Local Global Hybrid Order Optimal Tree G Optimal Dag <> Raw DP Global - Hybrid Order Figure 2.8: Test set log likelihood vs training time on Synthetic-15 network, d — 15, TV = 1500. The bottom figure presents the "good" algorithms in higher detail by removing the poor performers. Results for the factored model are an order of magnitude worse and therefore not plotted. Note that the DP algorithm actually took over two hours to compute. Adult (US Census) - 1 -0.08 -0.1 -0.12 J -* j ^ Gibbs Local / Global / Hybrid 1/ Y O Optimal Dag ( ) Raw DP 0 50 100 ISO 200 Time (seconds) Figure 2.9: Test set log likelihood vs training time on Adult dataset. d = 14, N = 49k. DP algorithm actually took over 350 hours to compute. The factored and maximum likelihood tree results are omitted since they are many orders of magnitude worse and ruin the graph's vertical scale. 21 Chapter 2. Structure learning by dynamic programming and MCMC Sachs (T-cell signalling pathway) I 7=*" 1 JL 3 -o.s • -0.6 • Global Hybrid Order <> Raw DP Gibbs Local Global Hybrid Order O Opimal Dag |' <^> Raw DP 150 Figure 2.10: Test set log likelihood vs training time on T-cel l dataset. d = 11, N = 5400. DP algorithm actually took over 90 hours to compute. The factored model and optimal tree plugins are again omitted for clarity. The bottom part of the figure is a zoom in of the best methods. 22 Chapter 2. Structure learning by dynamic programming and MCMC In addition to DP and the MCMC methods mentioned above, we also measured the performance of plug-in estimators consisting of: a fully factorized model (the disconnected graph), the maximum likelihood tree [7], and finally the MAP-optimal DAG gotten from the algorithm of [50]. We measure the likelihood of the test data as a function of training time, i(t). That is, to compute each term in I we use p(x,|D) = Yls=i P(xi\Cs), where St is the number of samples that can be computed in i seconds. Thus a method that mixes faster should produce better estimates. Note that, in the Dirichlet-multinomial case, we can quickly compute p(x\G") by plugging in the posterior mean parameters: p(xi^)=n^:r"G*=fc) ijk where Oijks = E[9ijk\D, G3}. If we have missing data, we can use standard Bayes net inference algorithms to compute p(x\Gs,6). In contrast, for DP, the "training" cost is computing the normalizing con-stant p(D), and the test time cost involves computing p(xi\D) = p{xi: D)/p(D) for each test case Xi separately. Hence we must run the DP algorithm m times to compute £ (each time computing the marginal likelihoods for all families on the augmented data set Xi,D). DP is thus similar to a non-parametric method in that it must keep around all the training data, and is expensive to apply at run-time. This method becomes even slower if x is missing components: sup-pose k binary features are missing, then we have to call the algorithm 2k times to compute p(x\D). For the first experiment, we generated several random networks, sampling the nodes' arities uniformly at random from between 2-4 and the parameters from a Dirichlet. Next, we sampled 100d records (where d is the number of nodes) and performed 10-fold cross-validation. Here, we just show results for a 15-node network, which is representative of the other synthetic cases. Fig-ure 2.8 plots the mean predictive likelihood across cross-validation folds and 5 independent sampler runs against training time. On the zoomed plot at the bottom, we can see that the hybrid and global MCMC methods are signifi-cantly better than order sampling. Furthermore, they seem to be better than 23 Chapter 2. Structure learning by dynamic programming and MCMC exact D P , which is perhaps being hurt by its modular prior. A l l of these Bayes model averaging ( B M A ) methods (except Gibbs) significantly beat the plugin estimators, including the MAP-op t ima l structure. In the next experiment we used the "Adult" US census dataset, which con-sists of 49,000 records with 14 attributes, such as "education", "age", etc. We use the discretized version of this data as previously used in [42]. The average arity of the variables is 7.7. The results are shown in Figure 2.9. The most accurate method is DP , since it does exact B M A (although using the modular prior), but it is also the slowest. Our D P + M C M C method (with B = 0.1) pro-vides a good approximation to this at a fraction of the cost (it took over 350 hours to compute the predictive likelihood using the D P algorithm). The other M H methods also do well, while Gibbs sampling does less well. The plug-in D A G is not as good as B M A , and the plug-in Chow-Liu tree and plug-in fac-tored model do so poorly on this dataset that their results are not shown (lest they distort the scale). (These results are averaged over 10 M C M C runs and over 10 cross-validation folds.) Finally, we applied the method to a biological data set [48] which consists of 11 protein concentration levels measured (using flow cytometry) under 6 different interventions, plus 3 unperturbed measurements. 600 measurements are taken in each condition yielding a total dataset of N = 5400 records. Sachs et al. discretized the data into 3 states, and we used this version of the data. We modified the marginal likelihood computations to take into account the interventional nature of the data as in [8]. The results are shown in Figure 2.10. Here we see that D P gives the best result, but takes 90 hours. The global, hybrid and order samplers all do almost as well at a fraction of the cost. The local proposal and Gibbs sampling perform about equally. A l l methods that perform B M A beat the optimal plugin. 24 Chapter 2. Structure learning by dynamic programming and MCMC Adult - Training set likelihood trajectories G i b b s Orde r - 3 . 7 - 3 . 8 - 3 . 9 1 1 - 3 . 9 0 50 100 150 200 50 100 150 2 0 0 0 50 100 150 200 x ! 0 s H V b r i d - 3 . 7 - 3 . 8 - 3 . 9 0 - 3 . 7 2 - 3 . 7 2 5 0 50 100 150 2 0 0 x l os O r d e r & G l o b a l 50 100 150 2 0 0 - 3 . 7 3 1 0 5 10 15 2 0 25 Figure 2.11: Samplers' training set log likelihood trace plots on Adult dataset. 4 traceplots of training set likelihood for each sampler, starting from randomly initialized values. The bottom-right figure combines runs from the order and global samplers and shows the behaviour of the chains in the first 25 seconds. 25 Chapter 2. Structure learning by dynamic programming and MCMC 2.5.4 Convergence diagnostics In Figure 2.11 we show a traceplot of the training set marginal likelihood of the different methods on the Adult dataset. (Other datasets give similar results.) We see that Gibbs is "sticky", that the local proposal explores a lot of poor configurations, but that both the global and order sampler do well. In the bottom right we zoom in on the plots to illustrate that the global sampler is lower variance and higher quality than the order sampler.- Although the difference does not seem that large, the other results in this section suggest that the DP proposal does in fact outperform the order sampler. 2.6 Summary and future work We have proposed a simple method for improving the convergence speed of M C M C samplers in the space of DAG models. Alternatively, our method may be seen as a way of overcoming some of the limitations of the DP algorithm of Koivisto and Sood [31, 32]. The logical next step is to attempt to scale the method beyond its current limit of 22 nodes, imposed by the exponential time and space complexity of the underlying DP algorithm. One way forward might be to sample partitions (layers) of the variables in a.similar fashion to [37], but using our DP-based sampler rather than Gibbs sampling to explore the resulting partitioned spaces. Not only has the DP-based sampler been demonstrated to outperform Gibbs, but it is able to exploit layering very efficiently. In particular, if there are d nodes, but the largest layer only has size m, then the DP algorithm only takes 0(d2m) time. Using this trick, [32] was able to use DP to compute exact edge feature posteriors for d = 100 nodes (using a manual partition). In future work, we will try to simultaneously sample partitions and graphs given partitions. This is a non-trivial task because the DP algorithm marginalizes over structure. The method of [37], for example, requires DAG samples to estimate parameters associated with partitioning. 26 Chapter 3 Structure learning with uncertain interventions 3.1 Introduction The use of Bayesian networks to represent causal models has become increasingly popular [45, 51]. In particular, there is much interest in learning the structure of these models from data. Given observational data, it is only possible to iden-tify the structure up to Markov equivalence. For example, the three models X—>V—>Z, X*— Y<— Z, and X*— Y—>Z all encode the same conditional inde-pendency statement, X _L Z \ Y . To distinguish between such models, we need interventional (experimental) data [16]. Most previous work has focused on the case of "perfect" interventions, in which it is assumed that an intervention sets a single variable to a specific state (as in a randomized experiment). This is the basis of Pearl's "do-calculus" (as in the verb "to do") [45]. A perfect intervention essentially "cuts off" the influence of the parents to the intervened node, and can be modeled as a structural change by performing "graph surgery" (removing incoming edges from the intervened node). Although some real-world interventions can be modeled in this way (such as gene knockouts), most interventions are not so precise in their effects. One possible relaxation of this model is to assume that interventions are "stochastic", meaning that they induce a distribution over states rather than a specific state [33]. A further relaxation is to assume that the effect of an intervention does not render the node independent of its parents, but simply 27 Chapter 3. Structure learning with uncertain interventions changes the parameters of the local distribution; this has been called a "mech-anism change" [52, 53] or "parametric change" [17]. For many situations, this is a more realistic model than perfect interventions, since it is often impossible to force variables into specific states. Here, we propose a further relaxation of the notion of intervention, and consider the case where the targets of intervention are uncertain. This extension is motivated by problems in systems biology and drug target discovery, where the effects of various chemicals that are added are not precisely known. In particular, each chemical may affect a hidden variable, which can in turn affect multiple observed variables, often in unknown ways. We model this by adding the intervention nodes to the graph, and then performing structure learning in this extended, two-layered graph. Our contributions are four fold. First, we show how to combine models of intervention — perfect, imperfect and uncertain — with a recently proposed al-gorithm for efficiently determining the exact posterior probabilities of the edges in a graph [31, 32], Second, we show empirically that it is possible to infer the true causal graph structure, even when the targets of interventions are un-certain, provided the interventions are able to affect enough nodes. Third, we apply our exact methodology to T-cell data that had previously been analyzed using M C M C [19, 48] and show that our uncertain intervention model is the best density estimator. Fourth, we utilize uncertain interventions to identify gene targets of cancer on the childhood acute lymphoblastic leukemia (ALL) data gathered by [58] and analyzed in [12, 58]. We believe our method is the first well-principled application of Bayesian networks to drug/disease target dis-covery. 3.2 Models of intervention We will first describe our probability model under the assumption that there are no interventions. Then we will describe ways to model the many kinds of interventions that have been proposed in the literature, culminating in our 28 Chapter 3. Structure learning with uncertain interventions (a) (b) Figure 3.1: Intervention model: None, (a) Plate notation, (b) the same model unrolled across 4 data cases. model of uncertain interventions. This will serve to situate our model in the context of previous work. 3.2.1 No interventions For the intervention-free case, we will assume that the conditional probabil-ity distribution (CPD) of each node in the graph is given by p(Xi\Xai,d,G) = fi(Xi\Xot, di), where Gi are the parents of i in G, &i are i's parameters, and /*() is some probability density function (e.g., multinomial or linear Gaussian). For the parameter prior p(9\G), we will make the usual assumptions of global and local independence, and parameter modularity (see [27] for details). We will fur-ther assume that each p(6i) is conjugate to /;, which allows for closed form com-putation of the marginal likelihood p{XUN\G) = /p{Xl:N\G,6)p(9)d0, where N is the number of data cases. For example, for multinomial-Dirichlet, the marginal likelihood for a family (a node and its parents) is given by [27] r N J n = l _ T T r ( « i j ) yr r ( q i j f c + NiJk) W r K + A ^ H r(« y f c) where Nijk = Z ! n = i Hxi = k,xG. = j) are the counts, and Nij = iZkNijk-(/(e) is the indicator function in which /(e) = 1 if event e is true and /(e) = 0 otherwise.) Also, ctijk are the pseudo counts (Dirichlet hyper parameters), 29 Chapter 3. Structure learning with uncertain interventions l Oi (a) (b) Figure 3.2: Intervention model: Perfect, (a) Plate notation, (b) the same model unrolled across 4 data cases. ctij = J2k, aijk, ri is the number of discrete states for Xt, and qt is the number of states for Xct- We will usually use the BDeu prior a,jk = l/<7i?"i [27]. (An anal-ogous formula can be derived for the normal-Gamma case [23].) The marginal likelihood of all the nodes is then given by p(X1:N\G) = YlLi P(X}:N\XGf), where d is the number of nodes. Figure 3.1 shows the non-interventional case as a graphical model. 3.2.2 Perfect interventions If we perform a perfect intervention on node i in data case n, then we set X™ = x*, where x* is the desired "target state" for node i (assumed to be fixed and known). We modify the CPD for this case to be p{Xi\XGi, 6) = I(X{ =x*). We see that Xi is effectively "cut off" from its parents Xa,L • Figure 3.2.(a) shows the perfect intervention model in plate notation, while (b) illustrates the idea on a local family with 4 data points. Namely, for i fixed, Figure 3.2.(b) "unrolls" the plate notation across 4 data. We see that in data cases 2 and 4 (marked in red) the perfect intervention was performed (I, = 1), cutting off Xi from its parents and corresponding parameters. Although not shown, the probability function over Xi's states has been collapsed onto the target state x*. 30 Chapter 3. Structure learning with uncertain interventions a. (a) (b) Figure 3.3: Intervention model: Imperfect, (a) Plate notation, (b) the same model unrolled across 4 data cases. X™ is node i in case n, Xg. are its parents. / " acts like a switching variable: If I™ = 1 (representing an intervention), then the hyper-parameters. 3.2.3 Imperfect interventions A simple way to model interventions is to introduce intervention nodes, that act like "switching parents": if If = 1, then we have performed an intervention on node i in case n and we use a different set of parameters than if If = 0, when we use the "normal" parameters. Specifically, we set p(Xi\XGi, h = 0,0, G) = fi(Xi\XGt,0?) and p(Xi\XGi,Ii = 1,0,G) = f^Xa^O]). (Note that the assumption that the functional form / ; does not change is made without loss of generality, since 0i can encode within it the specific type of function.) T ian and Pearl [52, 53] refer to this as a "mechanism change": see Figure 3.3. A special case of this is a perfect intervention, in which p(Xi\Xci, Ii = 1,0, G) = I{Xi = x*). To simplify notation, we assume every node has its own intervention node; if a node i is not intervenable, we simply clamp If = 0 for all n. When we have interventional data, we modify the local marginal likelihood formula by partitioning the data into those cases in which Xi was passively Xi uses the parameters 0\\ If / " = 0, then Xi uses the parameters 0®. o/i at are 31 Chapter 3. Structure learning with uncertain interventions Figure 3.4: Intervention model: Imperfect with unreliable extension. Compare to Figure 3.3. We can optionally add another switch node i?™, which can be used to model the degree of effectiveness of the intervention. observed, and those in which Xi was set by intervention: f] p(x?\xGi,8°)}p(6°)d6» /! l=0 n p(x?\xCi,o})]P(oi)d0i In the case of perfect interventions, this second factor evaluates to 1, so we can simply drop cases in which node i was set by intervention from the computation of the marginal likelihood of that node [8]. 3.2.4 Unreliable interventions An orthogonal issue to whether the intervention is perfect or imperfect is the reliability of the intervention, i.e., how often does the intervention succeed? One way to model this is to assume that each attempted intervention succeeds with probability <j>i and fails with probability 1 — cAjj this is what Korb et al. [33] call the degree of "effectiveness" of the intervention. We can associate a latent binary variable Rf to represent whether or not the intervention succeedeed or failed in case n, resulting in the mixture model P(Xi\xGi,ii = i,e,G) = X>(ifc = r)p{Xi\Xc^h = l,Ri = r,6,G) r = faMXiiXo^e}) + (l - &)/i(Xi|x G„0?)- (3-i) 32 Chapter 3. Structure learning with uncertain interventions X%. Xl i!=oi /?=ii /?=oi /f=ii - A ' ; Uv? '--X? HXf Wl (a) (b) Figure 3.5: Intervention model: Soft, (a) Plate notation, (b) the same model unrolled across 4 data cases. Proposed by [38]. x* is the known state into which we wish to force node i when we perform an intervention on it; Wi is the strength of this intervention. a\, the hyper-parameters of 9}, are a deterministic function of a1/, x* and Wi. Figure 3.4 illustrates the idea on the unreliable intervention model. An unreli-able, but otherwise perfect, intervention is modeled by setting p(Xi\Xat,Ii = l,Ri = 1,0,G) = I(Xi = x'i). Unfortunately, computing the exact marginal likelihood of a data case now becomes exponential in the number of R variables, because we have to sum over all 2' f il latent assignments. Although Figure 3.4 adds the indicator Rf to the imperfect model only, any of the other models of intervention under discussion could be augmented with the unreliable assumption also. 3.2.5 Soft interventions Another way to model imperfect interventions is as "soft" interventions, in which an intervention just increases the likelihood that a node enters its target state x*. Markowetz et al. [38] suggest using the same model of p(X i |A 'c7 i , 9, G) as before, but now the parameters 9° and 9} have dependent hyper-parameters. In particular, for the multinomial-Dirichlet case, O^j1 ~ Dir{a"'f), they assume the deterministic relation a^. = a?•. +Wiit, where j indexes states (conditioning cases) of XGt, t = x* is the target value for node i, et = (0,. . . ,0 ,1 ,0 , . . . ,0) 33 Chapter 3. Structure learning with uncertain interventions with a 1 in the t'th position, and u>i is the strength of the intervention. As Wi—>oo, this becomes a perfect intervention, while if Wi = 0 it reduces to an imperfect intervention. If the intervention strength wt is unknown, Markowetz et al. suggest putting a mixture model on wi} but it may be more appropriate to use the Ri mixture model mentioned above, where an intervention can succeed or fail on a case by case basis. Figure 3.5 shows the model graphically using plate notation. 3.2.6 Uncertain interventions Finally we come to our proposed model for representing interventions with un-certain targets, as well as uncertain effects. We no longer assume a one to one correspondence between intervention nodes /; and "regular" nodes Xi. In-stead, we assume that each intervention node /; may have multiple regular children. (Such interventions are sometimes said to be due to a "fat hand", which "touches" many variables at once.) If a regular node has multiple inter-vention parents, we create a new parameter vector for each possible combination of intervention parents: see Figure 3.6 for an example. We are interested in learning the connections from the intervention nodes to the regular nodes, as well as between the regular nodes. We do not allow connections between the intervention nodes, or from the regular nodes back to the intervention nodes, since we assume the intervention nodes are exogenous and fixed. We enforce these constraints by using a two layered graph structure, V = XUX, where X are the regular nodes and X are the intervention nodes. The addition of X motivates new notation, since the augmented adjacency matrix has a special block structure. The full adjacency matrix, denoted by H , is comprised of the intervention block F containing X nodes, and the backbone block G comprised of X nodes: We call the elements of F "target edges" since they correspond to edges X—>X 34 Chapter 3. Structure learning with uncertain interventions and the elements of G "backbone edges". As we will see in Section 3.3.1, the block structure of H reduces the time complexity of the DP algorithm that we use to perform exact Bayesian inference. To explain how we modify the marginal likelihood function, we need some more notation. Let Xai be the regular parents of node i, and Id be the inter-vention parents. Let 9\ be the parameters for node i given that its intervention parents have state L Then the marginal likelihood for a family becomes It is crucial that we assume that the interventions have local (albeit un-known) effects, otherwise they would not help us resolve Markov equivalency. To see this, note that if the distribution after an intervention, call it p(X|f?1), is unrelated to the distribution before an intervention, p(X\9a), then the overall marginal likelihood becomes a product of standard marginal likelihoods (gotten by integrating out 9° and 91). This gives us more data, but does not help us learn the causal structure (see [53] for more details). To see this, let us consider a simpler scenario (inspired by the analysis of [53]) in which we have a single intervention node Ii with target ch(Ii) = I. Suppose we observe iVo cases in which 7; = 0 and N\ cases in which Ii = 1. Let Nfjk be the counts in the first batch, N^k be the counts in the second batch, and Nijk = N°jk + N^k. Let Gx be the graph induced by the regular nodes. If the post interventional distribution is unconstrained (i.e., the parameters that generated the second batch of data are unrelated to the first set of parameters), then we get p ( ^ N \ x ^ J d N ) n P(X1:Nl,XNl+1:N2\ll-Nl = 0, J f 1 + 1 : / V 2 = 1,GX) rKQ -rrn*ijk + N?jk) • i1 ^ r ( a l j + ^ ) i i r(aiJk) 35 Chapter 3. Structure learning with uncertain interventions Figure 3.6: Example of a "Fat hand" intervention. Intervention 1 affects nodes 2 and 3, intervention 2 affects node 3. The parameters for node 3 are 0%]2{kj), where h = i, h =j, X2 = k and X3 = £. which is just the regular BDe likelihood applied to a larger dataset. But if we constrain the intervention to only affect node £, we get The unconditional marginal likelihood is then given by the mixture distribution In Section 3.3 we present an efficient way to compute this mixture, even in the case where there are multiple intervention nodes, each with potentially multiple targets. 3.2.7 The power of interventions The ability to recover the true causal structure (assuming no latent variables) using perfect and imperfect interventions has already been demonstrated both p{Xi]Nl,XNl + hN2\ll:N' =0,1-r(qijk + N i j k ) r(ao',fc) n<*i,i,k + Nijik) p(X\I,Gx) = Y,p(ch(Ii) = e)p(X\I, Gx,£)-e. 36 Chapter 3. Structure learning with uncertain interventions A A / \ / B C B C \ S \ w \ A A ; / \ / \ I B C B C | \ / \ _ \ / \ f Figure 3.7: Markov equivalence classes on the Cancer network, (a) The Cancer network, from [22]. (a-d) are Markov equivalent, (c-g) are equivalent under an intervention on B. (h) is the unique member under an intervention on A. Based on [53]. theoretically [16, 17, 52, 53] and empirically [8, 39, 52, 53, 57]. Specifically, each intervention determines the direction of the edges between the intervened nodes and its neighbors; this in turn may result in the direction of other edges being "compelled" [4]. For example, in Figure 3.7, we see that there are 4 graphs that are Markov equivalent to the true structure; given observational data alone, this is all we can infer. However, given enough interventions (perfect or imperfect) on B, we can eliminate the fourth graph (d), since it has the wrong parents for B. Given enough interventions on A, we can uniquely identify the graph, since we can identify the arcs out of A by intervention, the arcs into D since it is a v-structure, and the C-*E arc since it is compelled. In general, given a set of interventions and observational data, we can identify a graph up to intervention equivalence (see [52] for a precise definition). In Section 3.4.1, we will experimentally study the question of whether one can still learn the true structure from uncertain interventions (i.e., when the targets of intervention are a priori unknown), and if so, how much more data one needs compared to the case where the intervention targets are known. 37 Chapter 3. Structure learning with uncertain interventions 3.3 Algorithms for structure learning The Bayesian approach to structure learning avoids many of the conceptual problems that arise when trying to combine the results of potentially inconsis-tent conditional independency tests performed on different ("mutated") models [15]. In addition, it is particularly appropriate when the sample sizes are small, but "soft" prior knowledge is available, as in many systems biology experiments. However, for large structure learning problems Bayesian inference becomes in-tractable, and we will still have to resort to approximate and/or point estimation methods. 3.3.1 Exact algorithms for p(Hjj) and H*MAP On problems of less than d ~ 22 variables, we use the algorithm of Koivisto and Sood [31, 32], which can compute the exact posterior probabilities of all edges, p(Hij), using dynamic programming in 0(d2d) time, as we discussed in Section 2.2. The inputs to this algorithm are a "prior" over node orderings q-i(Ui), a "prior" over possible parent sets, pj(Gj), and a local marginal likelihood function for every node and every possible parent set, p(Xi\XGi)- They were also described in Chapter 2. Though we increase the effective number of nodes in using the uncertain intervention model, the block structure of H permits a reduction in the time complexity of the DP algorithm. Let dj = | J | be the number of intervention nodes, and dx = 1^ 1 be the number of regular nodes. The time complexity of the DP algorithm in this case is 0{d2dx + dk+1C{N)), where d = di + dX) and C(N) is the cost of computing each local marginal likelihood term. Note that layering is crucial for efficiently handling uncertain interventions, otherwise the algorithm would take 0(d2d) instead of 0{d2dx) time. Silander and Myllmaki recently devised an elegantly simple algorithm to compute the globally optimal DAG [50], which takes o(d22d~2) time and expo-nential space. The premise is that the optimal DAG must have a sink (a node without outgoing edges), which, by optimality, must have parents that score 38 Chapter 3. Structure learning with uncertain interventions the highest amongst all possible sets thereof. If this best sink and its incoming edges are fixed and removed from further consideration, the remaining nodes and edges must be optimal in the same fashion. When carried out recursively, these steps will construct the globally optimal DAG. We use this algorithm to compute the exact MAP, H*MAP. 3.3.2 Local search for H M / I P For much larger problems, Bayesian model averaging becomes infeasible to carry out, even by approximation methods. For the 271 variable A L L data (see Sec-tion 3.4.3), we will approximate the M A P by finding a structure ~H.MAP by local search that fits the training data well. As with BMA, there also exist two flavours of local search, which differ by the space they operate on: structure or order. In structure search, given a starting DAG, the algorithm considers all neighbouring DAGs that differ by a single edge addition, deletion or rever-sal then greedily chooses the one that most increases the training set marginal likelihood [27]. Recently, Schmidt et al. [49] compared the two approaches on several datasets and found that DAG search consistently outperformed order search on the structure learning task. Neighbour pruning is an optional pre-processing step that tests for independencies between variables, and rules out many structures or orders a priori. [49] found that if an appropriate neighbour pruning method is applied, DAG search scores better on training and test set likelihood as well. Therefore, we adopt the structure-space local search, forcing it to restart whenever a local maximum is reached, and taking only the highest scoring DAG based on training set marginal likelihood. 3.3.3 Iterative algorithm for VLMAP The block-H notation introduced in Section 3.2.6 is suggestive of a coordinate ascent algorithm, which learns an UMAP by alternating between learning F given G fixed, then G given F fixed. Even if the number of nodes dx is large, with a target edge fan-in constraint it will still be feasible to determine FMAP-39 Chapter 3. Structure learning with uncertain interventions This algorithm would seem particularly attractive for applications where the goal is to learn F only (G may be irrelevant). We implemented such an al-gorithm, and tested it according to the same procedures used in Section 3.4.1. However, we found that it becomes consistently trapped in local maxima, espe-cially when the globally optimal DAG algorithm of [50] is used for both steps. The culprit is in fixing F; when we do so, we effectively break the interventional nature of the data and create new local maxima. If the algorithm was allowed to look a sufficient number of steps into the future, it could escape these maxima; however, the extra computation required for the look-ahead would defeat the algorithm's original purpose. The problem can be somewhat ameliorated by using stochastic local search rather than finding the MAP, though in practise it appears to be preferable to simply learn F and G jointly. 3.4 Experimental results We first present some results on synthetic data generated from a Bayesian net-work of known structure, and then present results on a real biological dataset. 3.4.1 Synthetic data In this section, we experimentally study the question of whether one can still learn the true structure, even when the targets of intervention are a priori unknown, and if so, how much more data one needs compared to the case where the intervention targets are known5. We assessed this using the following experimental protocol. We considered the Cancer network shown in Figures 2.2 and 3.7, and then generated random multinomial CPDs by sampling from a Dirichlet distribution with hyper-parameters chosen by the method described in [6] (outlined in Section 2.5.2). For simplicity, we used binary nodes. We then 5 Tian and Pearl [52] briefly mention the case of "unknown focal variables" (which we are calling uncertain targets of intervention) in the context of constraint based learning methods, but do not present any algorithms for identifying focal variables. We are not aware of any other papers that address this question. 40 Chapter 3. Structure learning with uncertain interventions generated data using forwards sampling; the first 2000 cases D0 were from the original model, the second 2000 cases D\ from a "mutated" model, in which we performed a perfect intervention either on A or 5, forcing it to the "off" state in each case. Next we tried to learn back the structure using varying sample sizes of N £ {100,500,2000}. Specifically we used N observational samples and N interventional samples, D = (DQ:N, D\:N). We ran the algorithm using data D and under increasingly vague prior knowledge: (1) using the perfect inter-ventions model; (2) using the soft interventions model6; (3) using the imperfect model; and (4) using the uncertain interventions model. In the latter case, we also learned the children of the intervention node. As a control, we also tried just using observational data, D = DQ2N. Our results for the perfect and uncertain models are shown in Figure 3.8. On this network, the imperfect and soft intervention models perform very similar to the perfect case, though they require more data to achieve the same result. We see that with observational data alone, we are only able to recover the v-structure B—*D<—C, with the directions of the other arcs being uncertain (e.g., P(C—*E) w 0.75.) With perfect interventions on B, we can additionally recover the A-+B arc, and with perfect interventions on A, we can recover the graph uniquely, consistent with the theoretical results in Section 3.2.7. With uncertain interventions, we see that the entropy of the posterior on the regular edges is higher than when using perfect interventions, but it too reduces with sample size. Eventually the posterior converges to a delta function on the intervention equivalence class. We obtain similar results with other experiments on random graphs. This suggests that our proposed mechanism is able to learn causal structure even from uncertain interventions. Next, we turned our attention to the larger synthetic network "Cars Diag-nosis" introduced by [26], shown in Figure 3.9. Here, we will contrast the struc-tural recovery abilities of the perfect and uncertain models Receiver Operating 6 [38] do not discuss how to set the pushing strength Wi. We set it equal to 0.5N, so that the data does not overwhelm the hyper-parameter ctL f c-41 Chapter 3. Structure learning with uncertain interventions •|" Ground Truth N = 20 O A B C D E H = 6.31 N = 50 H = 5.40 N = 500 N = 2000 H = 1.86 H = 1.49 A B C D E H = 4.11 H = 2.53 H = 0.45 H = 0.09 0.2 0.4 0.6 0.8 1 Figure 3.8: Perfect vs uncertain interventions on the Cancer network. Results of structure learning on the Cancer network (Figure 3.7). Left column: ground truth. Subsequent columns: posterior edge probabilities p(Gy = 1\D) for increasing sample sizes N, where dark red denotes 1.0 and dark blue denotes 0.0. H is the entropy of the factored posterior T]._ p(Gij\D). See text for details. This figure is best viewed in colour. 42 Chapter 3. Structure learning with uncertain interventions Figure 3.9: Cars network ground truth, d = 12 model for car malfunction troubleshooting used by [26]. By selecting the appropriate two intervention nodes, marked here in red, it is possible to uniquely recover the structure. Characteristic (ROC) curve analysis, as in [31]. Again, we assign multinomial CPDs according to the method [6]. Without interventional data, it would be impossible to learn the orientation of some edges on the Cars Diagnosis network; however, if we intervene on nodes 1 and 9, it is possible to uniquely recover the original structure. Referring to Figure 3.9 we can see that under these two in-terventions, all edges become compelled that were not in the observational case. In this experiment, we account for sampling variability by creating 10 datasets for each of N = (20 ,200, 2000). For a given N, half of the data is observational and half interventional, with interventions generated according to the perfect model. Next, we attempted to learn the structure back for all 10 sets of data per sample size N, assuming both the perfect and uncertain models. Figure 3.10 presents the results using ROC plots to illustrate the effects of sampling different data sets. Since the DP algorithm outputs edge feature marginal probabilities, it is necessary to threshold them to induce hard decisions on edge presence. Let G? denote the marginal thresholded by 6 £ [0 ,1 ] , and G | be the corresponding ground truth. We say that G ^ is a True Positive {TP) if Gf- = 1 A G% = 1, or else if G^ = 1 A Gfj = 0, then Gei} is labeled a False Negative (FN). Similarly, G?* = 0 A G ^ = 0 means Gfj is a True Negative (TN), else if G^ = 0 A Gei}• = 1, we say that G^ is a False Positive (FP). The four classes {TP, FP,TN, FN} are mutually exclusive and exhaustive. Next, 43 Chapter 3. Structure learning with uncertain interventions let #TP indicate the sum of true positives over with identical definitions for the other classes. We now define Sensitivity = # T P / ( # T P + #FN), Specificity — #TN/(#FN+#TN) and InverseSpecificity = 1-Specificity. A ROC curve shows the inherent trade-off between specificity and sensitivity as the threshold 9 is varied. ROC plots are often summarized by a single quantity known as the AUC, which is simply the area under the curve. A perfect predictor receives an AUC of 1, while a random predictor gets 0.5. Given that the interventions are perfect in the generating network, we ex-pect the matching assumption for learning to perform best. As we can see from Figure 3.10, this is the case on average but the difference is not significant. As the sample size is increased, the gap between uncertain and perfect inter-vention performance closes. The AUC scores for the target edges suggest that learning the target edges requires less data, and happens before the backbone is determined. In Section 3.4.3, we analyze a d = 271 dataset for which Bayesian inference is intractable and instead use local search to approximate the MAP, yielding H-MAP- However, it is not clear how much is lost when estimating structure via approximate plugin (local search), versus exact plugin (MAP), versus exact B M A . Therefore, we repeated the synthetic Cars experiments, but this time comparing B M A , M A P and local search, all under the uncertain intervention model. For sample sizes of N 6 20,100,400,2000, we generated 20 independent datasets and ran each algorithm to obtain H as edge marginal probabilities (BMA) or single-DAG estimates (MAP and local search). Local search was allowed to run for 20 seconds (3 times as long as it took to compute the B M A or MAP) . Since ROC curves do not exist for hard assignments, we instead chose to show the specificity and sensitivity of each estimate. For B M A , we selected the equal-error rate point to threshold the marginals. This is the threshold for which the specificity equals the sensitivity, a fair measure of performance for this experiment. The results in Figure 3.11 suggest that except for extremely low sample sizes, point estimates are adequate estimates of both the backbone and target edges. Furthermore, local search performs within one standard deviation 44 Chapter 3. Structure learning with uncertain interventions Structure Recovery Perfect Interventions Structure Recovery Uncertain Interventions Mean AUC N = 20 N = 200 1 _ 0.6 I « I HI 1 N = 2000 0 Inv. specificity ' 0 Inv. specificity 1 JufcH 111 Figure 3.10: Structural recovery performance of perfect vs uncertain intervention models on Cars network. Each row of the figure denotes a different sample size (JV = 200 means there were 100 observational and 100 interventional data cases). The first two columns contain ROC curves showing the trade-off between specificity and sensitivity in thresholding the edge feature marginals produced by perfect and uncertain interventions. There are 10 curves per plot, corresponding to the 10 datasets generated per sample size. In the right column, the ROC curves have been summarized as AUC. The "Perfect" and "Uncertain" bars show structure recovery performance for the graph backbone (G), while "Target edges" applies only to the uncertain model, and is the AUC on intervention to backbone edges only (H). 45 Chapter 3. Structure learning with uncertain interventions of the exact M A P estimate. 3.4.2 T-cell data We now apply our methodology to a real biological dataset, which had pre-viously been analyzed using M C M C by Sachs et al. [48] (who used multiple restart simulated annealing in the space of DAGs), Werhli et al. [57] (who used Metropolis-Hastings in the space of node orderings), and Ellis and Wong [19] (who used equi-energy sampling in the space of node orderings). The purpose of our experiment is to determine the exact posterior over edges, and hence to assess the quality of the M C M C techniques, and also to learn the effects of the interventions that were performed. The dataset consists of 11 protein concentration levels measured under 6 dif-ferent interventions, plus 3 unperturbed measurements. The proteins in question constitute part of the signaling network of human T-cells, and therefore play a vital role in the immune system. See Figure 3.13(a) for a depiction of the commonly accepted "ground truth" network, including hidden nodes. The data in question were gathered using a technique called flow cytome-try, which can record phosphorylation levels of individual cells. This has two advantages compared to other measurement techniques: first, it avoids the in-formation loss commonly incurred by averaging over ensembles of cells; second, it creates relatively large sample sizes (we have N = 5400 data points in total, 600 per condition). The raw data was discretized into 3 states, representing low, medium and high activity. We obtained this discretized data directly from Sachs; see Fig-ure 3.12 for a visualization. This constituted the input to our algorithm. We tried two different analyses. In the first version, we assumed that the targets of intervention were known, and we modeled these using perfect inter-ventions (as did Sachs et al). The results are shown in Figure 3.13(c). These should be compared with the results of the M C M C analysis of Sachs et a l , which are shown in Figure 3.13(b), and the ground truth network, which is shown in 46 Chapter 3. Structure learning with uncertain interventions Performance on backbone Performance on target edges Specificity Sensitivity Specificity Sensitivity Figure 3.11: Structural recovery performance of exact B M A vs point estimation on Cars network. The leftmost two columns show performance on the graph backbone (G) while the rightmost two columns show it for the target edges (F). The results in each row were aggregated from 20 independent datasets with a sample size of TV (shown). TV = 200 means 100 interventional data and 100 observational data. The bottom row summarizes the plots above as curves. Note that since we chose the equal error rate threshold for B M A , we would expect the associated bar heights to be equal. This is not the case because there is no threshold setting such that Specificity(O) = Sensitivity(9) exactly. We note that the equal error rate threshold also corresponds to the maximum value taken on by f(9) = Specificity(9) + Sensitivity(9) therefore we take the maximum of f(9) as our operating point. 47 Chapter 3. Structure learning with uncertain interventions Training set Observed Biomolecule Figure 3.12: T-cell dataset. Discretized, 3-state data from [48]. Columns are the 11 measured proteins, rows are the 9 experimental conditions, 3 of which are "general stimulation" rather than specific interventions. The name of the chemical that was added in each case is shown on the right. The intended primary target is indicated by an E (for excitation) or I (for inhibition). This figure is best viewed in colour. -18 Chapter 3. Structure learning with uncertain interventions Figure 3.13(a). While there is substantial agreement between the three models, there are also many differences. For example, the ground truth shows no edge from jnk to p38, or from mekl2 to jnk, yet both inference methods detect such an edge. This may be due to the presence of various hidden variables. Looking at the data in Figure 3.12, mekl2 and jnk seem quite highly correlated, although this is obviously not enough evidence to suggest there should be an edge between them (as shown in [48], nearly all of the variables are significantly pairwise correlated). There are also several edges in our model that seem to be absent in the M C M C analysis of Sachs et al. This is possibly because Sachs et al. only perform model averaging over a "compendia of high scoring networks", as found by 500 restarts of simulated annealing, whereas our method averages over all graphs, and hence may detect support for many more edges. (Note that averaging over many sparse, but different, graphs can result in a dense set of marginal edge probabilities.) Also, the two methods use different graph priors p(G), and hence cannot be directly compared. We also applied the algorithm of Silander and Myllmaki [50] to compute the M A P DAG. On this dataset, it coincides exactly with Figure 3.13(c), suggest-ing that the posterior is highly peaked around the M A P structure. Since the algorithm does not sum over orderings, the bias discussed in Section 2.3 is not incurred, and the output M A P estimate is with respect to a uniform prior over DAGs. In the second experiment, we added the intervention nodes to the graph and learned their children, rather than pre-specifying them. The results are shown in Figure 3.13(d). We successfully identified the known targets of all but one of the 6 interventions. (We missed the G06967 —> pkc edge.) However, we also found that the interventions have multiple children, even though they were designed to target specific proteins. Upon further investigation, we found that each intervention typically affected a node and some of its immediate neighbors. For example, from the ground truth network in Figure 3.13(a), we see that Psitect 49 Chapter 3. Structure learning with uncertain interventions Figure 3.13: T-cell dataset results. Models of the biological data, (a) A partial model of the T-cell pathway, as currently accepted by biologists. The small round circles with numbers represent various interventions (green = ac-tivators, red = inhibitors). From [48]. Reprinted with permission from A A A S . (b) Edges with marginal probability above 0.5 as estimated by [48], (c) Edges with marginal probability above 0.5 as estimated by us, assuming known perfect interventions, (d) Edges with marginal probability above 0.5 as estimated by us, assuming uncertain, imperfect interventions, and a fan-in bound of fc = 2 for the target edges. The intervention nodes are in red, and edges from the intervention nodes are light gray. This figure is best viewed in colour. 50 Chapter 3. Structure learning with uncertain interventions • • l i •v*SK iiiptfil i : Perfea '• Soi Imperfect' Interventio n model (a) (b) Figure 3.14: Predictive likelihood scores of intervention models on T -cell dataset. Negative log predictive likelihood. Lower is better, (a) M A P es-timate (b) Exact B M A . Result obtained across 10-fold validation. Note: though difficult to see on this scale, exact B M A outperforms the plug-in estimate for uncertain interventions. (designated 8 in that figure) is.known to inhibit pip2; in our learned network (Figure 3.13(d)), we see that Psitect connects to pip2, but also to plcy, which is a neighbor of pip2. This is biologically plausible, since some of these interventions actually work by altering hidden variables, which can therefore cause changes in several neighboring visible variables. Also, although we missed the G06967 —> pkc edge, the other children of G06967 (plcy, pka, mekl2, erk and p38) seem to be strongly affected by G06967 when looking at the data in Figure 3.12. We also computed the M A P DAG and again found it to be identical to the thresholded edge marginals gotten by DP. We also tried analysing the continuous data using linear-Gaussian Bayes nets [23]. Following [19], we took a log transform of each variable and then standardized them. Our results are similar to [19], but our graph is much denser, suggesting that their M C M C scheme failed to visit sufficiently many modes. (Although once again our results are not directly comparable due to the different prior.) The graphs inferred using the Gaussian and multinomial models have much in common, but they also differ in many of the details. It is difficult to rigorously assess the quality of the obtained graphs when 51 Chapter 3. Structure learning with uncertain interventions there is no ground truth. (The biological model in Figure 3.13(a) is unlikely to be the "true" model that generated the data in Figure 3.12. Also, it contains hidden variables, so is not directly comparable to what we are learning.) The approach taken by Ellis et al. [19] was to compare the predictive log-likelihood in a cross-validation framework. This can also be done using the DP algorithm, by computing p(x\D) = p(x, D)/p(D); these normalization constants can be obtained by running the "forwards" algorithm of [32] using the "dummy" feature / = 1. Using 10-fold cross-validation we carried this procedure out for all of the intervention models. We also computed the predictive log-likelihood for the M A P structure under each model. Given the M A P DAG, this amounts to determining its posterior mean parameters, then evaluating the likelihood of each test point in the fold. Our results are shown in Figure 3.14. We see that our uncertain intervention model is the clear victor in both cases, suggesting that the assumption of perfect interventions may be poor for the T-cell data. Running time on T-cell data Experiments were performed on a laptop with a 2 GHz Intel Core Duo Processor and 2GB R A M running under Windows XP. For the 3-state T-cell data, with d = 11 nodes (using perfect interventions and a fan-in constraint of k = 5) and N = 5400, our Matlab implementation took 3.6 seconds to compute the marginal likelihood terms, while the DP algorithm took 0.4 seconds and the M A P DAG algorithm needed 0.8 seconds. For the case where we learned the effects of interventions (so d = 17), it took about 3.6 minutes (using a fan-in bound of = 5 for backbone edges and fc = 2 for target edges) to obtain the likelihood terms, 15 seconds for B M A and 80 seconds for the MAP. By comparison, the multiple restart simulated annealing approach used by Sachs et al. took several days. 52 Chapter 3. Structure learning with uncertain interventions 3.4.3 ALL data In this section we explore a promising application of the uncertain intervention model: drug/disease target discovery in gene networks. The complex interac-tions in most gene networks are poorly understood, but there is hope that they can be estimated from measurements such as gene expression microarray data. Learning a protein signalling network under normal cellular conditions was the goal in Section 3.4.2, but a more interesting, if more difficult, application is to determine the effect of external agents like drugs or disease on gene interactions. A biologist may or may not want to learn the backbone network as well. Since our intervention model directly enables this type of analysis, we apply it to the gene expression data obtained by Yeoh et al. [58], which was previously analyzed in [12, 58]. The data is comprised of measurements of 12,000 genes, from 327 humans suffering from different forms of acute lymphoblastic leukemia (ALL). A L L is a heterogeneous cancer, meaning that it is manifested by several subtypes that vary by their genetic influence and consequently in their response to treatment. The dataset of [58] contains 7 classes, 6 of which represent the common A L L subtypes HYPERDIP > 50, E2A-PBX1, B C R - A B L , T E L - A M L 1 , M L L and T - A L L , and the final class aggregating several less common subtypes. We followed [12] and omitted all but 271 genes from our analysis, using the Chi-square-based filtering method of [58] which selects the top 40 discriminative genes for each subtype (9 genes were chosen multiply chosen across subtypes, yielding 271 unique genes). We also discretized the data to 3 levels "under-expressed" (+1), "unchanged" (0) and "overexpressed" (—1) using the same procedure as [12]. Specifically, if (/ij,fjj) were the mean and standard deviation of gene i, values less than pi — <7j were mapped to —1, greater than p^ + a* to + 1 and the remainder to 0. The resulting dataset is shown in Figure 3.15.(a). Dejori et al. assume that each A L L subtype acts on a single gene, which then causes other genes downstream in the network to change. They learned a Bayesian network using simulated annealing search on the 327-case dataset, ignoring the fact that the data samples come from different conditions. Let 53 Chapter 3. Structure learning with uncertain interventions Sampled data (a) (b) Figure 3.15: A L L dataset. (a) Training dataset (b) data sampled from model under same conditions as training set (cancer subtype orders and frequencies). This figure is best viewed in colour. Dk denote all data cases representing ALL subtype index k ranging over the 7 subtypes. Using exact inference they generated a sample D[ = {x ~ p ( X _ j | X ; = +1)} and another = {x ~ p(X-i\Xi = —1)} for each gene i, and compared the Euclidean distance between D\ and Dk for each k. They declared the cause of type k cancer to be the gene that, when overexpressed, generated data D[ that most closely resembled Dk- Note that they use a Bayesian network simply as a density estimator for discrete data; they make no attempt to interpret the structure of the graph. When they set Xi = +1 or X-i = —1, they treat this as an observation, rather than a Pearl-style do-action, so they could have in principle used any other kind of density estimator. There are two obvious shortcomings to the approach of [12]. First, if we believe that the presence of cancer gives rise to changes in the gene interactions, then it is not sensible to learn a single Bayesian network across multiple cancer conditions. The analysis of Dejori et al. avoids this issue by not discussing the graph structure they learn. However, they do so implicitly when they use their Bayesian network to simulate data from "no cancer" condition. Secondly, from a computational standpoint, their method is very expensive due to the need for 54 Chapter 3. Structure learning with uncertain interventions Highest scoring DAG Target edges (a) (b) Figure 3.16: Highest scoring D A G found on A L L dataset. Highest scoring found across 20 parallel, 24 hour-long runs of restarting local search, (a) DAG backbone, G , (b) corresponding target edges, F . Next-best structure had a score approximately e 6 4 times lower. inference on a large network to sample from p(X_j |Xj = ± 1 ) . This problem would be greatly compounded if they considered the cause of A L L to be the mutation of b < B genes, since this would require inference to be performed 0((d)) times per subtype. In our approach, we augment the 271 backbone variables with 7 binary intervention nodes encoding the presence or absence of the A L L subtypes and attempt to learn H . Following [12], we assume that each gene has at most one cancer subtype parent. The problem size is well beyond the limitations of the exact DP or M A P algorithms, therefore we use local search, modified to support uncertain interventions. We tried the M M P C neighbour pruning algorithm introduced in [56] as a way improve search, but found that even though unrestricted DAG search is slower by a factor of 4, it finds better-scoring graphs. M M P C likely gave poor results because it identifies sets of potential neighbours (parents and children) using conditional independency tests that are undoubtably unreliable with a sample size of 371 for d = 271 variables. We 55 Chapter 3. Structure learning with uncertain interventions No Cancer Figure 3.17: Predicted gene-expression profile for a patient without ALL. The procedure of [12] is arguably incapable of estimating these profiles. Their Figure 3.(b) appears to support this claim. ran 20 searches seeded at random initializations for 24 hours each, and allowed them to restart at a new random DAG whenever they became caught in a local maximum. In the results that follow, we adopt the highest scoring DAG H^g as our plugin estimate of the structure, H . Figure 3.16 shows the graph. We note that Hxs shared substantial similarity with other high scoring graphs. Since our model is generative we can sample data from it and then com-pare this data to the training data to determine if the model is sensible. Fig-ure 3.15.(b) displays data sampled from R*£,s under the same conditions as the training set, meaning that in the cases where, for example, HYPERDIP > 50 was present in the training data, we clamped the subtype's correspond-ing intervention node to "on" for the same cases in the synthetic dataset. It is apparent that our model has captured much of the detail from the original data. In Figure 3.17 we sample 327 cases from the non-interventional condition. This condition.was not contained in the training data, but our model is capable of learning it through the assumed locality of interventions. Figure 3.18.(a)-(b) shows the expression profile for 327 simulated patients with subtype E2A-PBX1 or M L L . Our model can also estimate expression profiles for hypothetical pa-tients who have unluckily contracted more than one subtype of A L L . One such 56 Chapter 3. Structure learning with uncertain interventions SO 100 ISO 200 2S0 300 Data casts (c) Figure 3.18: Sampled gene-expression profile for a patient wi th A L L subtype (a) E 2 A - P B X 1 or (b) M L L or (c) both. 57 Chapter 3. Structure learning with uncertain interventions combination, E2A-PBX1 and M L L , is shown in Figure 3.18.(c), which appears to be reasonable after referring to (a) and (b). The main objective of this analysis is to determine the genetic targets of A L L . These are easy enough to read off of F L S and can be seen in Figure 3.16.(b). As can be seen, the model has picked up many targets for each subtype, especially for T E L - A M L 1 . We require a method to rank these edges in order to compare our results with [12], who report the top 5 scoring gene targets per A L L subtype. One method would be to analyze the top M scoring - DAGs found by local search, and assign a score to all of the target edges in their union according to the number of times each edge appeared. If the DAGs were samples from the posterior, this would be valid; however, greedy local search is certainly not sampling from the posterior. Instead, we adopt a more principled approach by computing, for each target edge in FLs, the cost of removing that edge with all other structure being fixed. The expression for this weight, denoted W^x with i g l a n intervention node and x e X & backbone gene node, is given by: Wi,x = \ogp(D\Gx,i^x) - \ogp(D\Gx), (3.2) where Gx G X are the backbone parents of x and p{D\-) are local marginal likelihood terms. WjtX is the plugin estimate to a Bayes factor that answers the same question, but integrates over the remaining structure rather than holding it fixed. We can also compute WiiX on the edges not chosen by FLS, with the quantity now representing the gain in adding these edges. By the fact that local search find a local maximum, WilX is certain to be positive for edges turned "on" in FLS and negative otherwise. With Wi/X we can plot a spectrum of target edge scores across the 271 genes, for each subtype. These are shown in Figures 3.19-3.21.(a). Note that these represent "monogenic activations" only. We could also compute W~i:X for vector arguments of i, but we follow [12] and only report results on single gene perturbations. Also, many of the negative WiiX are not shown as their score is negative infinity. These infinite values arise from the initial assumption that genes can only be affected by one cancer. Therefore, sparsity/density on the plot 58 Chapter 3. Structure learning with uncertain interventions shows which genes do not have any parent in the intervention (cancer subtype) layer. These figures also show where our results overlap with [12]; genes which Dejori et al. chose as their top 5 are shown in red with a black square instead of a circle. Next to the squares is a number between 1 and 5 that indicates how highly [12] ranked the gene to be a cause of A L L . Referring to Figure 3.19, we see that our top-scored gene for subtypes E2A-PBX1 and B C R - A B L agree with Dejori et al.'s. This is a positive result, because these genes are proto-oncogenes suspected to cause those A L L subtypes. Our respective spectra for subtype E2A-PBX1 also closely agree, though it is difficult to take the comparison any further, since we do not have access to their detailed results. Another interesting capability of our model is in determining the type of interaction between the cancer and the genes it targets. Given F/,s we compute its posterior mean parameters and then examine the conditional probability table for a particular gene target Xi that has a cancer parent Ii. We marginalize the backbone parents, and look at the parameters corresponding to the "on" state of the cancer: p ( X l | / l = 1 ) = E x „ x G i p ( ^ | x G i , / t = i)-If most of the mass resides in the "overexpressed" state, p(X, = = 1), we say the cancer is excitory for that gene, while if the "underexpressed" state, p(Xi = — 1|J, = 1), dominated we would say it is inhibitory. We plot the ex-pression of target genes in Figures 3.19-3.21.(b). Red upwards arrows indicate excitation, while blue downwards arrows show inhibition. The remaining prob-ability mass corresponding to the "no change" state is not shown, but can be easily derived from the fact that the three probabilities sum to unity. 3.5 Summary and future work We have shown how to apply the dynamic programming algorithm of Koivisto and Sood [31, 32] to learn causal structure from interventional data. We then introduced the model of uncertain interventions, which enables the discovery of 59 Chapter 3. Structure learning with uncertain interventions E2A's effect on directly connected genes (b) BCR's effect on directly connected genes 25 271 (c) (d) Figure 3.19: A L L subtypes E2A-PBX1 and B C R - A B L results. E2A-PBX1: (a)-(b), BCR-ABL: (c)-(d). Left: Target edge strength WIA for chosen edges (positive) and absent edges (negative). Points in red marked by a square indicate overlap with the results of [12]. Gaps correspond to target edges which are impossible. The result for E2A-PBX1 agrees strongly with the corresponding spectrum of [12]. Right: Cancer's effect on gene's expression level. This analysis is sensible only for target edges that were chosen (i.e. present in FLS)-GO Chapter 3. Structure learning with uncertain interventions Cancertype; HYPERDIP>50 HYPER DIP >50's effect on directly connected genes tl t 11 1 ' p yi |F * 25 271 (c) Cancertype: MLL ! ? T {> 1 6 2 (b) TEL's effect on directly connected genes (d) MLL's effec! on directly ccnnected genes (e) Cancer type; T-ALL (0 T-all's effect on directly connected genes 180 25 271 (g) ( h ) Figure 3.20: A L L subtypes Hyperdip > 50, TEL-AML1, M L L and T - A L L results. Hyperdip > 50: (a)-(b), T E L - A M L 1 : (c)-(d), M L L : (e)-(f), T -ALL: (g)-(h). Left: Target edge strength Wi]X Right; cancer's effect on expression level. ( i l Chapter 3. Structure learning with uncertain interventions Cancef type: OTHERS OTHERS's effect on directly connected genes | HI U— -0.5 • ~~~ \ 45 90 135 180 225 271 1 4 5 90 135 1 80 225 271 (a) (b) Figure 3.21: Other A L L subtypes' results, (a) Target edge strength WilX (b) cancer's effect on expression level. No corresponding result to (a) is presented in [12], presumably because this class does not represent a homogeneous A L L subtype. the target of an intervention. Leveraging the dynamic programming algorithm, we established that little data is needed to learn intervention targets, while only modestly more is required to simultaneously determine the graph backbone. We applied this model to two gene expression datasets that are of great relevance to modern systems biology, and demonstrated novel results in both cases. In the near future, we will apply the uncertain intervention model to the data of Ideker et al. [30], which consists of gene expression data from the galactose pathway of the yeast Saccharomyces cerevisiae. Using another methodology, Hallen et al. [25] report learning the gene targets of the compound Galactose; we expect to reproduce this result using our approach. Another interesting extension would be to apply the uncertain intervention idea to the active learning case [44, 54], where one has to decide which inter-ventions to perform. „? 9 1 I I ! 1 A w 62 Bibliography [1] A. -L. Barabsi and Z. N . Oltvai. Network biology: Understanding the cell's functional organization. Nature Reviews Genetics, 5:101-113, 2004. [2] G. Brightwell and P. Winkler. Computing linear extensions is #P-complete. In STOC, 1991. [3] W. Buntine. Theory refinemement on Bayesian networks. In UAI, 1991. [4] D. Chickering. A transformational characterization of equivalent Bayesian network structures. In UAI, 1995. [5] D. Chickering, D. Heckerman, and C. Meek. A Bayesian Approach to Learning Bayesian Networks with Local Structure. In UAI, 1997. [6] D. Chickering and C. Meek. Finding Optimal Bayesian Networks. In UAI, 2002. [7] C. K. Chow and C. N . Liu. Approximating discrete probability distributions with dependence trees. IEEE Trans, on Info. Theory, 14:462-67, 1968. [8] G. Cooper and C. Yoo. Causal discovery from a mixture of experimental and observational data. In UAI, 1999. [9] R. G. Cowell, A. P. Dawid, S. L. Lauritzen, and D. J. Spiegelhalter. Prob-abilistic Networks and Expert Systems. Springer, 1999. [10] D. Dash and G. Cooper. Model Averaging for Prediction with Discrete Bayesian Networks. J. of Machine Learning Research, 5:1177-1203, 2004. 63 Bibliography [11] N. de Freitas, P. Hjen-Srensen, M . I. Jordan, and S. Russell. Variational M C M C . In UAI, 2001. [12] Mathaus Dejori and Martin Stetter. Identifying interventional and pathogenic mechanisms by generative inverse modeling of gene expression profiles. Journal of Computational Biology, 11(6):1135-1148, 2004. [13] D. Eaton and K. Murphy. Bayesian structure learning using dynamic pro-gramming and M C M C . In UAI, 2007. [14] D. Eaton and K. Murphy. Exact Bayesian structure learning from uncertain interventions. In A I/Statistics, 2007. [15] F. Eberhardt. Sufficient condition for pooling data from different distri-butions. In First Symposium on Philosophy, History, and Methodology of Error, 2006. [16] F. Eberhardt, C. Glymour, and R. Schemes. On the number of experiments sufficient and in the worst case necessary to identify all causal relations among N variables. In UAI, 2005. [17] F. Eberhardt, C. Glymour, and R. Scheines. Interventions and causal in-ference. In 20th Mtg. Philos. of Sci. Assoc., 2006. [18] D. Edwards. Introduction to graphical modelling. Springer, 2000. 2nd edition. [19] B. Ellis and W. Wong. Sampling Bayesian Networks quickly. In Interface, 2006. [20] J. Friedman, T. Hastie, and R. Tibshirani. Additive logistic regression: a statistical view of boosting. Annals of statistics, 28(2):337-374, 2000. [21] N. Friedman and D. Roller. Being Bayesian about Network Structure: A Bayesian Approach to Structure Discovery in Bayesian Networks. Machine Learning, 50:95-126, 2003. 64 Bibliography [22] N . Friedman, K. Murphy, and S. Russell. Learning the structure of dynamic probabilistic networks. In UAI, 1998. [23] D. Geiger and D. Heckerman. Parameter priors for directed acyclic graph-ical models and the characterization of several probability distributions. The Annals of Statistics, 30(5):1412-1440, 2002. [24] P. Giudici and R. Castelo. Improving Markov chain Monte Carlo model search for data mining. Machine Learning, 50(1—2):127—158, January 2003. [25] K. Hallen, J. Bjorkegren, and J. Tegner. Detection of compound mode of action by computational integration of whole genome measurements and genetic perturbations. BMC Bioinformatics, 7(51), 2006. [26] D. Heckerman, J. Breese, and K. Rommelse. Troubleshooting under uncer-tainty. Technical Report MSR-TR-94-07, Microsoft Research, 1994. [27] D. Heckerman, D. Geiger, and M . Chickering. Learning Bayesian net-works: the combination of knowledge and statistical data. Machine Learn-ing, 20(3):197-243, 1995. [28] D. Husmeier. Sensitivity and specificity of inferring genetic regulatory in-teractions from microarray experiments with dynamic Bayesian networks. Bioinformatics, 19:2271-2282, 2003. [29] K. -B . Hwang and B.-T. Zhang. Bayesian model averaging of Bayesian net-work classifiers over multiple node-orders: application to sparse datasets. IEEE Trans, on Systems, Man and Cybernetics, 35(6):1302-1310, 2005. [30] T. Ideker, V. Thorsson, J. Ranish, R. Christmas, J. Buhler, R. Bumgarner, R. Aebersold, and L. Hood. Integrated genomic and proteomic analysis of a systematically perturned metabolic network. Science, 2001. Submitted. [31] M . Koivisto. Advances in exact Bayesian structure discovery in Bayesian networks. In UAI, 2006. 65 Bibliography [32] M . Koivisto and K. Sood. Exact Bayesian structure discovery in Bayesian networks. J. of Machine Learning Research, 5:549-573, 2004. [33] K. Korb, L. Hope, A. Nicholson, and K. Axnick. Varieties of causal inter-vention. In Pacific Rim Conference on AI, 2004. [34] D. Madigan, J. Gavrin, and A. Raftery. Enhancing the predictive per-formance of Bayesian graphical models. Communications in Statistics -Theory and Methods, 24:2271-2292, 1995. [35] D. Madigan and A. Raftery. Model selection and accounting for model uncertainty in graphical models using Occam's window. J. of the Am. Stat. Assoc., 89:1535-1546, 1994. [36] D. Madigan and J. York. Bayesian graphical models for discrete data. Intl. Statistical Review, 63:215-232, 1995. [37] V . Mansinghka, C. Kemp, J. Tenenbaum, and T. Griffiths. Structured priors for structure learning. In UAI, 2006. [38] F. Markowetz, S. Grossmann, and R.Spang. Probabilistic soft interventions in Conditional Gaussian networks. In 10th Al/Stats, 2005. [39] F. Markowetz and R. Spang. Evaluating the effect of perturbations in reconstructing network topologies. In Proc. 3rd Intl. Wk. on Distrib. Stat. Computing, 2003. [40] M . Meila and T. Jaakkola. Tractable Bayesian learning of tree belief net-works. Statistics and Computing, 16:77-92, 2006. [41] M . Meila and M . I. Jordan. Learning with mixtures of trees. J. of Machine Learning Research, 1:1-48, 2000. [42] Andrew Moore and Weng-Keen Wong. Optimal reinsertion: A new search operator for accelerated and more accurate Bayesian network structure learning. In Intl. Conf. on Machine Learning, pages 552-559, 2003. 66 Bibliography [43] Andrew W. Moore and Mary S. Lee. Cached sufficient statistics for efficient machine learning with large datasets. J. of AI Research, 8:67-91, 1998. [44] K. Murphy. Active learning of causal Bayes net structure. Technical report, Comp. Sci. Div., UC Berkeley, 2001. [45] J. Pearl. Causality: Models, Reasoning and Inference. Cambridge Univ. Press, 2000. [46] C. Robert and G. Casella. Monte Carlo Statisical Methods. Springer, 2004. 2nd edition. [47] R. W. Robinson. Counting labeled acyclic digraphs. In F. Harary, editor, New Directions in the Theory of Graphs, pages 239-273. Academic Press, 1973. [48] K. Sachs, O. Perez, D. Pe'er, D. Lauffenburger, and G. Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308, 2005. [49] M . Schmidt, G. Fung, and R. Rosales. Generalized smooth L l regulariza-tion. In Intl. Conf. on Machine Learning, 2007. Submitted. [50] T. Silander and P. Myllmaki. A simple approach for finding the globally optimal Bayesian network structure. In UAI, 2006. [51] P. Spirtes, C. Glymour, and R. Schemes. Causation, Prediction, and Search. MIT Press, 2000. 2nd edition. [52] J. Tian and J. Pearl. Causal discovery from changes. In UAI, 2001. [53] J. Tian and J. Pearl. Causal discovery from changes: a Bayesian approach. Technical report, UCLA, 2001. [54] S. Tong and D. Roller. Active learning for structure in Bayesian networks. In Intl. Joint Conf. on AI, 2001. 67 Bibliography [55] I. Tsamardinos, L. Brown, and C. Aliferis. The max-min hill-climbing Bayesian network structure learning algorithm. Machine learning, 2006. [56] I Tsamardinos, L Brown, and C Aliferis. The max-min hill-climbing Bayesian network structure learning algorithm. Machine Learning, 65(l):31-78, 2006. [57] A. Werhli, M. Grzegorczyk, and D. Husmeier. Comparative evalua-tion of reverse engineering gene regulatory networks with relevance net-works, graphical Gaussian models and Bayesian networks. Bioinformatics, 22(20):2523-2531, 2006. [58] EJ Yeoh, M E Rossa, SA Shurtleff, and W K Williams et al. Classifica-tion, subtype discovery, and prediction of outcome in pediatric acute lym-phoblastic leukemia by gene expression profiling. Cancer Cell, 1:133-143, 2002. 68 Part II Appendices Appendix A Software The code developed and used for this thesis is freely available at w w w . c s . u b c . c a / - m u r p h y k / S t r u c t u r e L e a r n i n g , and is accompanied by documentation for the major features. The package is named Bayesian Network Structure Learning (BNSL); it was written in Matlab, and consequently can be run on any platform supported by that software. BNSL has no other external dependencies. BNSL is capable of all the exact and approximate Bayesian structure learn-ing methods reported in this thesis. Both static and dynamic networks can be learned, using multinomial or linear-Gaussian CPDs. Al l four models of inter-vention can be used, along with purely observational data. Table A . l enumerates the major functionality Name Type Limitation Dynamic programming [31] Exact d « 20 - 22 Exhaustive enumeration Exact d = 6 Gibbs sampling [37] Approx. d » 20 Order space sampling [37] Approx. dc* 20 Structure space M C M C [13, 24, 35] Approx. Depends (*) Optimal DAG [50] Exact d 2 20 - 22 Local search [27] Approx. d = 500 Table A . l : Major features of BNSL software (*) If the global proposal is used, the limitations are the same as dynamic programming; however, if pure local moves are made, the technical limitation becomes d ~ 50 (in practice, local moves will not mix well for d > 10). 70
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use.
Graduate
Bayesian network structure learning for the uncertain experimentalist : with applications to network biology
Text