It is now widely understood that most diseases have a genetic component and yet do not follow the simple Mendelian inheritance patterns of dominant or recessive traits [1]. Presumably, these diseases result from the interacting effects of multiple genes in combination with one or more environmental risk factors [2]. Such complex phenotypic traits are characterized by a high level of unpredictability, as both the number and nature of interactions are difficult to distinguish using conventional methods. Nevertheless, uncovering the genetic basis for disease and deciphering the relative contribution of environmental exposure are critical steps toward the goal of developing an effective system of “personalized” medicine [3, 4].

Machine learning methods, in particular Bayesian networks (BNs), have the potential to help disentangle the web of relations among genes, environment, and disease [5]. BNs are a multivariate modelling method able to simultaneously account for gene-gene (epistasis) and gene-environment interactions, as well as leverage the diagnostic potential of clinical or physiological factors [6]. BNs also lead directly to prognostic models: a constructed network, however complex, can be used to efficiently compute the probability that an individual with a particular genotype and environmental exposure will exhibit the phenotype of interest. BNs have been applied in a variety of settings for the purposes of causal study and probabilistic prediction, including medical diagnosis, crime and terrorism risk, forensic science, and ecological conservation (see [7]). In bioinformatics, they have been used to analyze gene expression data [8, 9], derive protein signaling networks [10–12], predict protein-protein interactions [13], perform pedigree analysis [14], conduct genetic epidemiological studies [5], and assess the performance of microsatellite markers on cancer recurrence [15]. In this paper, we review the potential for BNs to contribute to revealing the genetic and environmental basis of disease.

where x = {x ₁ ,…x _n } are the variables (nodes in the BN) and θ = {θ ₁ ,…,θ _n } are the BN’s parameters, where each θ _i is the set of parameters necessary to specify the distribution of the variable x _i given its parents pa(x _i ). (We adopt the notation of [18] in which lower case symbols are used to indicate particular realizations of the corresponding uppercase variables.) Importantly, the factorization given by eq. (1) allows complex systems to be analyzed and modeled using a limited number of local relationships as building blocks.

If each of the five variables in Figure 1 is assumed to be binary, this factorization reduces the number of parameters (conditional probabilities) required to specify the full joint distribution from 64 to 12. This allows large networks to be parameterized from relatively small data sets while still capturing the kinds of interactions that are typical of complex traits.

A number of recent reviews [20–22] emphasize the complex nature of the gene-environment-disease relationship and call for new statistical methods to address such complexities. The parametric linear modeling framework (e.g. logistic regression) traditionally used in association studies, while familiar and easy to implement and interpret, is limited in its ability to detect interactions [23]. Studies employing this approach also typically start by considering only one SNP at a time and then include interaction effects only for those SNPs that exhibit independent marginal effects [22]. This procedure ignores the broader environmental context and multigenic nature of complex disease. It also runs the risk of identifying too many associated SNPs as a result of linkage disequilibrium or evolutionarily induced dependencies.

These concerns are exemplified by the relationships shown in Figure 1. A simple one-at-a-time search procedure might identify associations between C, E, and S1, as long as the environmental exposure and SNP both have direct effects on cancer. But if S1 serves primarily to mediate the effects of E on C, its role would be missed if interactions are not properly considered. Further, because of its association with S1, it is possible that S2 would be identified as being associated with C. This would introduce a redundancy if S2 and S1 are in linkage disequilibrium, or a false association if S2 and S1 are not on the same gene or region of linkage disequilibrium [19]. Finally, traditional approaches focused exclusively on direct relations would likely overlook the roles of S3 and B, although their identification could have utility for disease prediction, identification, and prevention.

Improved procedures for implementing linear models can partially address these shortcomings by putting more emphasis on interactions [24] or efficient parameterizations [25]. Yet, the high numbers of higher-order interactions that are now believed to underlie most human disease are not compatible with the strengths of traditional statistical methods [26]. Multifactor dimensionality reduction (MDR) is one example of a novel computational strategy for detecting and characterizing multiple non-linear interactions in the absence of detectable marginal effects [27, 28]. MDR works by pooling genotypes from multiple loci to create new variables before subsequent associative analysis. A user-friendly open source software package for implementing MDR is freely available from http://www.epistasis.org.

Even after pre-processing by MDR, linear logistic regression-type methods require large sample sizes to achieve significant power. They also require focus on a single well-defined phenotype. In contrast, through efficient factorization, BNs can break down the discovery process of a complex system into separate investigations of smaller components [29]. This means disease phenotypes can be defined by multiple factors, including clinical or physiological diagnostic variables. Further, as opposed to creating problems of multicollinearity, the associations between candidate predictor variables are naturally accounted for when defining a BN’s conditional probability distributions [19]. In the next section, we describe the means by which BNs can be used in genetic association studies.

Learning a BN can refer to data-based inference of either: (i) the conditional probability parameters for a given structure or (ii) the underlying graphical structure itself. We focus here on the process of structural learning, as the discovery of novel relations between genes, the environment, and disease is typically considered a harder problem of greater interest than the estimation of effect strength for a known factor. Often, parameter values are estimated concurrently as part of structural learning.

When the data are complete, the modularity of BNs provided by the Markov condition (eq. 1) also facilitates the process of structural learning. This is because the likelihood function used as the basis for most scoring metrics can be decomposed into a product of terms, each of which depends only on a variable’s parents and the conditional probability parameters. This means that candidate networks generated by small changes (adding or reversing an edge) can be evaluated locally, without regard to changes made elsewhere in the network. This allows for efficient learning algorithms.

Unfortunately, when data are incomplete – as inevitably occurs with large SNP arrays and hence gene-environment-disease data – decomposition of the likelihood function is no longer possible. This makes structural learning in the presence of missing values a computational challenge. Most state-of-the-art algorithms, including the ones reviewed in the previous section, are not able to learn BNs from incomplete data. One solution might be to simply discard observations with missing values or replace them with the mean of the observations. However, in either case, statistical power is lost and the joint distribution may be distorted. Therefore, in this section we review the use of the popular expectation-maximization (EM) algorithm [43] for handling missing values in BN learning. We start with the situation of parameter learning for a known BN structure, followed by the more difficult problem of structural learning.

Thus far, we have implicitly assumed that the BNs being discussed obey causal semantics; that is, every directed edge represents an association between parent and child that can be assumed to be causal, as in Figure 1. This requires that the Markov condition of Eq. 1 is assumed to extend to the notion of causality; given the effects on a variable of its immediate causes, that variable is independent of all earlier causes. Of course, in the case of genetic data, this assumption may be indirect, as we do not necessarily suppose that SNPs themselves are causal to the disease but rather that a SNP that is in linkage disequilibrium with the measured SNP causes a functional alteration [21].

There are two more assumptions that are required in order for a BN to be interpreted as representing causality: (i) if any two variables are correlated, then it must be assumed that one is the causal parent of the other or there is a third variable causing both, and (ii) there must be no variables missing from the network that are causal parents of two or more of the variables contained in the network [18]. Therefore when any important variables of a system are left unobserved – as is often the case – it is not clear that a BN learned from data (i.e., found to be consistent with the distributions of observed values) will have a directed structure that can interpreted causally. In fact, for any single directed graph selected as “best” according to a score-based algorithm, there will generally be many other graphs that are consistent with same conditional independencies implied by the selected graph. As with the partially directed graphs emerging from constraint-based algorithms, these graphs are termed observationally equivalent and the set of graphs in the same equivalence class can be readily identified using existing algorithms [52].

Claiming that a BN represents causality is tempting because not only does a causal BN provide mechanistic insight into a system but it allows the effect of interventions (e.g., reductions in exposure, drug treatments) to be correctly predicted. Unfortunately, as stated above, the possibility of unobserved variables (often referred to as hidden or latent variables) and the existence of a typically large equivalence class usually preclude causal interpretation. On the other hand, interpreting BNs in terms of causality is not necessary for extracting meaningful information from learned structures. Schlosberg et al. [53] learn a fully directed BN involving a target phenotype and 533 SNPs using a hill climbing algorithm. Perhaps counter-intuitively, they interpret all descendants of the phenotype to be the causes, and compare the results to what would have been obtained had they interpreted the causes to be: (i) all members of the phenotype’s Markov blanket or (ii) only the children of the phenotype.

Sebastiani and Perls [19] provide the rationalization for the approach of Schlosberg et al. [53]: being a complex trait, one would expect the presence of disease to be modulated by a large number of SNPs. These SNPs may have modifying effects on disease status, meaning that the association between one SNP and the disease will affect the strength of association of another SNP. The consequence is that only a limited number of SNPs can be practically detected as parents of the disease in causal models that are statistically learned from data. Sebastiani and Perls suggest that this limitation is removed by the use of diagnostic models, in which we allow SNPs and environmental exposures to be children, rather than exclusively parents, of the disease. In this way, the ability to detect the association of each SNP with the disease is not influenced by the association of other SNPs. This type of structure also more accurately represents the data-generating mechanism of a case–control study in which subjects are chosen based on their disease status, rather than at random. In practice, because the other parents of a child of the disease also contain information on disease status, Sebastiani and Perls consider the entire Markov blanket of the disease node as potential causes. This is also the approach taken by Han et al. [54] using a constraint-based approach to structure learning. We will refer to such network representations as non-causal.

The probability of disease given a particular environmental exposure and genetic profile is not directly encoded in a non-causal BN, but it can be readily computed using Bayes’ theorem. This technique is used in the analyses of Sebastiani et al. [6] and Ramoni et al. [55]. Because of the counter-intuitive nature of such diagnostic models, Sebastiani and Perls [18] suggest presenting results as an undirected graph representing mutual associations, rather than causal dependencies. The undirected counterpart of a directed graph, referred to as a Markov network, is formed by connecting nodes with a common child and then removing the direction of all edges. This is equivalent to connecting each node to its Markov blanket [56].

To exemplify the various BN structural learning algorithms and causal interpretations, we attempt to learn the relationships between a limited number of gene-environment-disease variables collected as part of a population-based case–control study of bladder cancer in New Hampshire. In addition to the presence/absence of bladder cancer, the full dataset includes over 1477 SNPs in cancer-related genes, detailed smoking assessment, gender, age, possible environmental risk factors including arsenic exposure, and selected biomarkers [57, 58]. Informed consent was obtained from each participant and all procedures and study materials were approved by the Committee for the Protection of Human Subjects at Dartmouth College. In the present context, we use the BN learning method to further assess reported interactions between DNA repair genes and arsenic exposure in increasing bladder cancer risk [59]. Specifically, in assessing polymorphisms in the XRCC3 and ERCC2/XPD genes using logistic regression, evidence of an increased risk of bladder cancer among those in the top arsenic exposure decile was observed for those with a variant allele of the double-strand break repair gene XRCC3. Therefore, we focus our analysis on 11 variables that allow us to explore the role of SNPs in the gene XRCC3 at positions 03, 04, and 241 and ERCC2/XPD at positions 03, 09, and 312. Arsenic exposure is represented by toenail arsenic level as an internal biomarker, which, following Andrew et al. (2009), is dichotomized at the 90th percentile. We also include the known risk factors: gender, age (≤60 or >60), and smoking status (dichotomized as never and former/current). While this set of 11 variables is certainly very small for a gene association study, we believe this allows us to more clearly demonstrate the BN method.

For analysis, we use the open-source software package bnlearn[60] in the statistical and graphical environment R[61] (http://www.r-project.org). The bnlearn package implements all of the algorithms described in Subsection Common algorithmic implementations in addition to a number of their variants. It also supports parameter estimation for given networks, conditional probability queries, model comparison and manipulation, random data generation, and plotting.

Complete data

We begin by comparing the Grow-Shrink(GS), Incremental Association Markov Blanket (IAMB), Hill-Climbing (HC), and Max-Min Hill-Climbing (MMHC) algorithms to a subsample of our dataset consisting of 424 controls and 226 cases without any missing values over the variables of interest.

For the constraint-based GS and IAMB algorithms and the hybrid MMHC algorithm, we used the Pearson’s χ ² test as the basis for the conditional independence tests with an alpha value (nominal type I error rate) of 0.05. Other available options in bnlearn include a chi-square test on the mutual information, a shrinkage estimator for the mutual information, and an experimental AIC-based independence test (see [60] for details).

For the score-based HC algorithm and the hybrid MMHC algorithm, we chose log(K2) as the scoring metric, as it is intermediate to the AIC and BIC in its penalization of complexity. Other frequentist scoring options for discrete variables in bnlearn include the multinomial log-likelihood, AIC, and BIC. To avoid convergence on a local maximum, we implemented 100 random restarts, each with 5 edge perturbations.

To first limit our results to fully prognostic models, we used the available ‘blacklist’ option in bnlearn to disallow any unreasonable casual relationships (e.g. SMOKING as a parent of AGE, non-genetic variables as parents of SNPs, CANCER as a parent of GENDER or SMOKER, SNPs as parents of AGE, TOENAIL_AS as a parent of AGE, GENDER, or SMOKER).

Results obtained from the four algorithms are broadly similar. In fact, the two strictly constraint-based algorithms and the two algorithms including a scoring component each returned identical structures (Figure 2). In all networks, SMOKER and GENDER are identified as parents of CANCER, and GENDER is additionally a parent of SMOKER. All networks show AGE to be a parent of TOENAIL_AS. While none of the networks show AGE as a parent of CANCER, as might be expected, this is not surprising in the present context given that controls were selected in a manner that minimized differences in age relative to cases [57].

None of the structures learned from the complete data cases shows any relation between the selected SNPs and other variables. The linkage disequilibrium groups found are the same across algorithms, with the exception of an edge between XRCC3_03 and XPD_312 in the network produced by the HC and MMHC algorithms. The directions of edges differ slightly between the strictly constraint-based algorithms and those with a scoring component, but as these are understood to only represent associations rather than genuine causal relationships this is not a serious concern.

When the networks learned by the algorithms are scored post-hoc using the various criteria (Table 3), it can be seen that those returned by HC and MMHC yield higher scores on all metrics except the BIC. Computationally, the constraint-based algorithms performed fewer tests than the ones with a score component, demonstrating their computational efficiency. As the computational complexity is polynomial in the number of tests, this could be a significant consideration for datasets with many more variables.

	GS	IAMB	HC	MMHC
Tests	259	366	750	921
Directed arcs	10		8
Log-likelihood	-4045.78		-4043.24
AIC	-4073.78		-4072.24
log(K2)	-4125.07		-4124.27
BIC	-4136.46		-4137.16

The highest value of each score is in bold.

We next consider a non-causal network (i.e., one in which edges do not necessarily follow causal direction) by removing from the ‘blacklist’ those edges for which CANCER is a parent and re-running the score-based HC algorithm. We choose the HC algorithm because of the results of Table 3 and because, having dropped the assumption of causality, we are interested in high scoring networks, rather than those that accurately represent causal direction. The result (Figure 3) shows three variables: GENDER, SMOKER and XRCC3_241 to have direct relations with CANCER, the first as a parent and the second two as children. Additionally, as a parent of XRCC3_241, XRCC3_04 is in the Markov blanket of CANCER. Therefore, all four of these variables should be considered further as candidate contributory causes of bladder cancer.

Missing data

Andrew et al. [59] found evidence of gene-environment interaction between XRCC3 241 and high toenail arsenic levels on bladder cancer risk. Thus far, our BN analysis using only the 424 controls and 226 cases with complete data has not been able to confirm this association. To investigate the possible influence of missing values, we next employed the EM algorithm to enable us to use the full dataset of 665 controls and 448 cases. (Most of the missing values occur in the variables representing SNPs.) To simplify the search process, we only consider the hypothesis that TOENAIL_AS should be included in the Markov blanket of CANCER. This could occur if this variable is either a parent of CANCER, a child of CANCER, or another parent of a child of CANCER (Figure 4). Thus, including the network representing the null hypothesis, we consider five possible structures. After application of the EM algorithm to each structure, the expected values of the log-likelihood, AIC, log(K2), and BIC scores were calculated as the basis for model comparison.

Results indicate that the network including TOENAIL_AS as a parent of XRCC3_241 provides the best fit to the data according to the first three scoring criteria (Table 4). Only the BIC score favors the null hypothesis of no association. The inclusion of TOENAIL_AS in the Markov blanket of CANCER is a result that is qualitatively consistent with the findings of Andrew et al. [59]. Additionally, the fact that TOENAIL_AS is a common parent of XRCC3_241, together with CANCER and XRCC3_04, indicates the presence of gene-environment interactions in determining cancer risk (Figure 5).

Added edge	Log-Likelihood	AIC	log(K2)	BIC
none	-3824.1	-3837.1	-3861.7	-3869.7
TOENAIL_AS → CANCER	-3824.6	-3839.6	-3865.7	-3877.2
CANCER → TOENAIL_AS	-3824.6	-3838.6	-3865.3	-3873.7
TOENAIL_AS → SMOKER	-3820.9	-3837.9	-3864.6	-3880.5
TOENAIL_AS → XRCC3_241	-3815.8	-3832.8	-3859.6	-3875.4

The highest scoring network according to each score is in bold.

To quantitatively explore these interactions, the risk of an individual having bladder cancer given gender, smoking status, toenail arsenic level, and XRCC3 04 and XRCC3 241 variation can be computed using Bayes’ theorem as applied to the final network in Figure 5. We used the cpquery function in bnlearn to perform the inference required to calculate the conditional probabilities of cancer. We then computed odds ratios relative to a female, non-smoker, with low toenail arsenic levels and the wildtype allele at both the 04 and 241 positions on the XRCC3 gene. The relative odds ratios of some of the 32 possible combinations (Table 5) show that, when the known risk factors gender and smoking are each present alone, they are associated with comparable increases in bladder cancer risk (see cases 6 and 7). For females without other observed risk factors, the variant allele at position 04 is associated with higher bladder cancer risk than the double wildtype form of this gene, while a double-variant form is associated with lower risk (compare cases 1, 2, 3, and 4). The association between cancer risk and arsenic exposure depends on genotype. For any combination of gender and smoking status, the risk of bladder cancer is elevated when toenail arsenic levels are high for those with a variant at position 241 and wildtype at position 04 (compare cases 3, 5, 9 and 11 and cases 10 and 13). However, for those with a wildtype at position 241 and variant at position 04, the cancer risk does not increase with arsenic level (compare cases 12 and 13). A double variant genotype is associated with lower risk of bladder cancer than single variant genotypes (compare cases 1, 3 and 4, and cases 8, 12, and 14). Such results show the potential for BNs to capture the complex role of multiple SNPs, as well as their interplay with exposure and background variables in determining disease susceptibility. Traditional logistic regression results corresponding to the relations discovered in the final BN confirm the statistical significance of the single factors, as well as the interaction between XRCC3_241 and XRCC3_04 (Table 6). The combination of high arsenic exposure and XRCC3_241 variation reveals an elevated odds ratio but an insignificant p-value of 0.14, suggesting that further study is necessary to confirm this association.

Case	Gender	Smoker	Toenail As	XRCC3_241	XRCC3_04	Odds ratio	# of subjects
1	female	no	low	variant	variant	0.68	13
2	female	no	low	wildtype	wildtype	1.0 (ref)	22
3	female	no	low	variant	wildtype	1.10	16
4	female	no	low	wildtype	variant	1.45	39
5	female	no	high	variant	wildtype	2.02	2
6	female	yes	low	wildtype	wildtype	2.22	27
7	male	no	low	wildtype	wildtype	2.23	20
8	male	yes	high	variant	variant	2.25	5
9	female	yes	low	variant	wildtype	2.45	24
10	male	yes	low	variant	wildtype	4.36	181
11	female	yes	high	variant	wildtype	4.48	9
12	male	yes	high	wildtype	variant	5.17	14
13	male	yes	low	wildtype	variant	5.75	95
14	male	yes	high	variant	wildtype	7.99	14

	Coefficient	Std. error	P-value	Odds ratio	95% CI
Intercept	-1.382	0.170	<0.00001	1	-
S	0.640	0.154	0.00003	1.90	(1.40, 2.56)
G	0.633	0.143	<0.00001	1.88	(1.42, 2.49)
X4	0.403	0.179	0.024	1.50	(1.05, 2.12)
A:X241	0.435	0.295	0.140	1.55	(0.87, 2.76)
X241:X4	-0.731	0.240	0.002	0.48	(0.30, 0.77)

While our example application is quite small in comparison to genome-wide association studies, it exemplifies some important points regarding BN structural learning. First, while the claim that a BN can represent causal relations is tempting, strictly enforcing a causal interpretation by disallowing nonsensical causal directions can limit the identification of important associations. For example, while the influence of gender and smoking on bladder cancer was clearly identified by our causal model (Figure 2), the association of cancer with the XRCC3_241 SNP was missed. This may not be surprising, given that cancer is effectively the independent variable in the process of selecting subjects for a case–control study. In the non-causal model (Figure 3), the XRCC3_241 SNP is identified as a child of CANCER (as is SMOKER), while GENDER continues to be a parent of both CANCER and SMOKER. This is consistent with the role of GENDER as a consideration in the process of selecting controls.

Interestingly, in the non-causal model, the XRCC3_04 SNP is revealed as being associated with CANCER by virtue of its membership in the Markov blanket. The V-structure that is formed at the XRCC3_241 node is a distinctive feature of BN modeling; such a structure implies that the two parents are marginally (i.e., unconditionally) independent, but become dependent when conditioned on the value of the child. A chi-squared test applied to the data confirms the fact that XRCC3_04 and CANCER are marginally independent (χ² = 0.85, p = 0.36, n = 1113), but become significantly positively associated (χ² = 7.75, p = 0.005, n = 558) conditional on XRCC3_241 being wildtype. A one-at-a-time search strategy (e.g. Andrew 2009) will typically miss such an association. This is an example of Simpson’s paradox, and appropriately capturing such situations in a BN allows for the accurate representation of complex relations.

In most real-world datasets, much information is lost when only complete observations are considered in statistical analysis. The EM algorithm provides a practical means for estimating model parameters without disregarding observations with missing values. In our example, this greatly increased our sample size and allowed for the discovery of toenail arsenic levels as a significant predictor of bladder cancer. TOENAIL_AS was not determined to be a parent of CANCER, but rather was found to be a parent of XRCC3_241, together with CANCER and XRCC3_04. This indicates the presence of gene-environment interactions, as reflected in the pattern of odds ratios calculated from this structure.

By focusing our attention on the role of TOENAIL_AS and its potential membership in the Markov blanket of CANCER, we were able to perform a comprehensive comparison of candidate structures after implementing the EM algorithm on each. Normally, one would be interested in comparing too many different structures to apply EM and scoring algorithms to them all. In such cases, the SEM algorithm or one of its variations would be necessary. Structural learning in the presence of missing data continues to be an active area of research.

We employed a score-based method when learning our non-causal structure and when considering TOENAIL_AS as a member of CANCER’s Markov blanket. While constraint-based methods can be more efficient, they can also be sensitive to failures in hypothesis tests. Further, if one has already given up causal claims, then a high-scoring structure is more important than a network that accurately represents causal direction. Of course, the choice of scoring metric is influential in determining which network is “best.” Because we view the process of BN structural learning as primarily an exploratory phase of data analysis that should be followed by rigorous replication studies [62], we do not see a need to rely on scoring criteria, such as the BIC, that penalize harshly for the number of edges.

Our example contained all discrete or easily discretized variables. While we believe this is generally representative of problems concerning the genetic dissection of disease, it may be that results are sensitive to the particular discretization of continuous variables. For example, based on the analysis of Andrew et al. [59], we divided toenail arsenic levels into “low” and “high” based on the 90th percentile of the observed values. There may be reasons why other thresholds would be more predictive or biologically relevant. When there is no clear basis for choosing the thresholds, other investigators have used the quartile boundaries [19]. Continuous variables could also be kept as such, in which case the conditional probability distributions in eq. (1) are represented by conditional density functions [35]. This does not necessarily present a problem for learning BN parameters for a given structure, as statistical methods – including the handling of missing data – are readily available. However, structural learning becomes significantly more difficult when variables are continuous, as the number and type of possible dependence relations and interactions becomes infinite. However, under some assumptions, such as linear relations and conditionally normal distributions, effective algorithms have been worked out [63]. Using kernel density estimators to model the conditional distribution, Hofmann and Tresp [64] were able to eliminate the reliance on normality. Further research is being conducted on this practically relevant topic [38].

We constrained ourselves here to discussing networks containing only variables for which there are some recorded data. Of course, in most cases, not all relevant aspects of a problem have been observed. Such hidden variables can present a problem for network structural learning, as omission of nodes effectively amounts to marginalization of the underlying joint distribution, potentially leading to complex dependencies among the remaining, observed variables. For example, in the gene-environment-disease context, data may be available on several environmental biomarkers and health outcomes, as well as a number of predisposing genetic or sociocultural factors. The relation between these is presumably mediated by the level of exposure to an environmental stressor, rendering many of the effects and predisposing factors conditionally independent. However, if the actual exposure level is not being measured, then all the observed variables will appear to be related to each other, likely in complex ways. By explicitly including a node representing a hidden variable in a network – even if there are no recorded data on that variable – the learned models are likely to be simpler and less prone to overfitting [65]. Fortunately, the structural EM algorithm can handle hidden variables analogously to how it handles missing values of otherwise observed quantities [65]. However, how to choose the number and placement of hidden variables in a BN remains an active area of research.

Despite being applied to Bayesian networks, the process of structural learning as we have described it this far has not been truly Bayesian in spirit. That is, we have not incorporated prior knowledge regarding possible structures or attempted to calculate posterior probabilities of model structures or features (e.g. the presence of a particular edge or the value of a particular parameter). As mentioned in Section Methodological approaches, the BIC score represents an approximation of the full marginal likelihood, which is comparable to the Bayesian posterior with a vague prior. The K2 score also represents a special case of the Bayesian posterior resulting from a uniform Dirichlet prior. A further generalization, still assuming a Dirichlet distribution but with prior knowledge able to be incorporated as an equivalent sample size was developed by Heckerman et al. [66] and is referred to as the Bayesian Dirichlet equivalent (BDe) score. Friedman [65] describes a version of the structural EM algorithm that directly addresses Bayesian model selection.

Even when the model selection process is Bayesian, in that it incorporates prior knowledge, typically only the single structure with the greatest posterior probability is maintained for further prediction and inference. However, when the amount of data (i.e. number of observations) is small relative to the size of the network (i.e. number of nodes or edges), there are likely to be many structures that conform with the data nearly equally well [67]. In such a situation – which is the norm when working with gene-environment-disease data with many SNPs relative to the sample size – the choice of a single “best” structure is largely arbitrary. Another data manifestation from the same population could have led to a very different final model. Zhang [68] has developed a fully Bayesian graphical method for large-scale association mapping, called BEAM3, which yields posterior probabilities of association. Yet, in most situations there are exponentially many structures that provide a “reasonable” representation of the observed data, making comprehensive investigation or communication of all “good” structures (i.e. those with a non-negligible posterior probability) impossible. For this reason, Friedman and Koller [67] propose a method for computing the Bayesian posterior of a particular structural feature, defined as the total posterior probability of all structures that contain the feature. Such features might include the presence of a particular edge between two nodes, the choice of a node’s parents, or the Markov blanket of a node. In some cases, robust assessment of such features may be more relevant to biological discovery than full articulation of (potentially fragile) network structures.