Sensitivity of Bayesian Network Structure Learning to Search Algorithms and Scoring Functions: A Comparative Study Across Three Healthcare Datasets

1. Introduction

Clinical decision making increasingly relies on data driven models, yet the models that perform best on held out accuracy are often the hardest to explain to a physician who must justify a diagnosis. Bayesian networks remain attractive for this reason: a learned directed acyclic graph, together with conditional probability tables, is a scientific object that can be read and checked against domain knowledge, unlike a black box classifier (Zhang et al. 2024).

Learning a Bayesian network from data is a search and score problem: a search strategy proposes candidate graphs and a scoring function ranks them (Chickering 1996). Hill Climbing and Tabu Search are the most widely used search strategies; we add Greedy Equivalence Search (GES) as a third, structurally distinct strategy searching over equivalence classes of graphs (Chickering 1996). BIC and BDeu are the most widely used scoring functions; we add K2, a long established Bayesian criterion also used as a baseline elsewhere in this literature (Khanteymoori et al. 2008). Prior work applying several algorithms to a breast cancer Bayesian network reported accuracy differences across algorithms (Sun et al. 2023), raising the question of how much of that effect originates in the search procedure versus the scoring criterion, and whether such a pattern holds beyond one dataset. Related work has shown several algorithms unstable with respect to variable column order (Kitson and Constantinou 2025), a different, arbitrary source of instability from the one we study: holding data and variable order fixed, how much does structure change with the algorithm and with the scoring function, and how stable is either comparison under resampling, across independent datasets.

We address this question on three widely used healthcare datasets: Breast Cancer Wisconsin Diagnostic (Cruz Ramirez et al. 2007; Zhang et al. 2024), Pima Indians Diabetes, and UCI Cleveland Heart Disease, chosen because they are compact, commonly used structure learning benchmarks spanning a range of dimensionality and predictor type. For each dataset we apply redundancy aware feature selection where predictor count exceeds our interpretability threshold, compare three search algorithms under BIC, repeat under BDeu and K2, and test robustness to feature selection and discretization choices, and to resampling, on the reference dataset. Our contribution is a quantitative, three dataset answer, with a compact, reproducible pipeline reusable on further datasets. Section 2 reviews related work; Section 3 details methodology; Section 4 reports results; Section 5 discusses them; Section 6 and Section 7 address limitations, future work, and conclusions.

2. Related Work

2.1. Bayesian Networks in Clinical Research

Bayesian networks have been applied to breast cancer data for diagnosis and prognosis. Diagnostic studies report accurate but algorithm dependent classifiers (Cruz Ramirez et al. 2007; Sun et al. 2023; Zhang et al. 2024). Prognosis oriented work combines Bayesian networks with survival models or hybrid architectures, consistently recovering plausible predictors such as tumor stage (Toffaha et al. 2025; Bani Mohammad and Ahmad 2025; Su et al. 2023; Iraji et al. 2020; Shahmirzalou et al. 2023; Piot et al. 2024). None of these studies tests whether a different search algorithm or scoring function would produce a materially different graph, nor whether findings on one dataset replicate on another, so their interpretive claims are implicitly conditioned on choices whose consequences are rarely examined.

2.2. Search Strategies and Search Spaces

Hill Climbing and Tabu Search perform almost identically to each other and to a genetic algorithm across benchmark networks (Ribeiro and Maciel 2023), consistent with our results but obtained on synthetic networks without varying the scoring function. Other heuristics include ant colony and genetic search over orderings (Wu et al. 2010), memetic anytime search (Lee et al. 2017), and swarm based search under BDe (Kareem and Okur 2020; He et al. 2024). A separate literature varies the search space itself: equivalence class search (Chickering 1996), hybrid skeleton and greedy orientation (Tsamardinos et al. 2006), constrained or restricted graph spaces (Acid and de Campos 2003; Acid et al. 2013), expert knowledge restrictions (de Campos and Castellano 2007; Abramovici et al. 2008; de Santana et al. 2007), and population based metaheuristics (Blanco et al. 2003; Sahin et al. 2004; Biba 2009). Hill Climbing and Tabu Search have also been shown sensitive to variable column order, with stabilized variants proposed (Kitson and Constantinou 2025), a different source of instability than the one we study.

2.3. Scoring Functions

BIC has repeatedly outperformed BDeu and related scores at recovering true structures on data from real networks (Liu et al. 2012; Liao et al. 2022), behaving more consistently across dataset sizes than BDeu, which is sensitive to its equivalent sample size (Chua and Ong 2017; Khanteymoori et al. 2008). BDeu can violate a basic regularity property, adding an uninformative parent even at zero conditional entropy (Suzuki 2016). K2, one of the earliest Bayesian scores, has been used as a baseline alongside BIC and AIC in classification oriented comparisons (Khanteymoori et al. 2008). Genomic scale applications, at far higher dimensionality than here, include local gene network recovery (Prestat et al. 2013), proteomic biomarker discovery (Wang et al. 2023), and gene environment disease modeling (Su et al. 2013; Beretta et al. 2017).

Overall, algorithm comparisons are typically run under one fixed score and one dataset, and score comparisons under one fixed search, rather than evaluated jointly, across multiple datasets, with an explicit robustness check. This is the gap the present study addresses.

3. Methodology

3.1. Bayesian Network Representation

A Bayesian network over variables $X=\{X_1,\dots ,X_n\}$ pairs a directed acyclic graph $G$ with conditional probability tables, factorizing the joint distribution as

$\begin{array}{r}P\left(X_1,\dots ,X_n\right)={\displaystyle \prod _{i=1}^n}P\left(X_i\mid {\mathrm{P}\mathrm{a}}_G\left(X_i\right)\right),\end{array}$where ${\mathrm{P}\mathrm{a}}_G\left(X_i\right)$ denotes the parents of $X_i$ in $G$. All variables are discrete; conditional probability tables are estimated by maximum likelihood once $G$ is fixed.

3.2. Datasets and Preprocessing

We use three datasets, all publicly available from the UCI Machine Learning Repository. Breast Cancer Wisconsin Diagnostic (Wolberg et al. 1993) (569 patients, 30 continuous morphology predictors, binary diagnosis). Pima Indians Diabetes (Kahn 1994) (768 patients, 8 continuous predictors, binary outcome); five predictors encode missing measurements as zero, a well documented issue for this dataset, addressed by replacing zeros with the column median over non-zero values. UCI Cleveland Heart Disease (Janosi et al. 1989) (303 patients, 13 predictors, binary outcome), of which 5 (age, resting blood pressure, cholesterol, maximum heart rate, ST depression) are continuous and 8 already categorical and left undiscretized. Continuous predictors are quantile discretized into $k$ equal frequency bins ($k=3$ default; Section 4.5 checks sensitivity), preferred over equal width bins since several predictors are right skewed.

Where candidate predictors exceed our interpretability threshold (10), we apply redundancy aware feature selection. Breast Cancer Wisconsin’s thirty predictors report ten morphological properties at three summary statistics each, several near collinear by construction; we rank predictors by mutual information with the label,

$\begin{array}{r}I\left(X_j;Y\right)=\sum _x\sum _{y}P\left(x,y\right)\mathrm{l}\mathrm{o}\mathrm{g}{\displaystyle \frac{P\left(x,y\right)}{P\left(x\right)P\left(y\right)}},\end{array}$then greedily select among the twenty highest ranked candidates by maximizing

$\begin{array}{r}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\left(f\right)=I\left(f;Y\right)-\lambda \cdot {\displaystyle \frac{1}{\left|S\right|}}\sum _{s\in S}|\rho \left(f,s\right)|,\end{array}$where $S$ is the already selected set and $\rho$ is Pearson correlation on the original continuous measurements, retaining eleven predictors. We set $\lambda =1.3$ from an internal sweep over $\lambda \in \{0,0.5,1.0,1.3,1.5,2.0,2.5,3.0\}$ (Section 4.5): family diversity among selected predictors rises from 5 at $\lambda =0$ to a plateau of 8 spanning $\lambda \in \left[1.3,1.5\right]$, then falls as redundant variants are reintroduced; we use $1.3$, the plateau’s smaller boundary, reporting the sweep directly rather than an external default (Han et al. 2026). The same criterion applied to Heart Disease’s thirteen predictors retains nine; Pima’s eight are used without reduction.

3.3. Search Strategies

Hill Climbing and Tabu Search explore the neighborhood of single edge additions, removals, and reversals that preserve acyclicity, starting from an empty graph with at most four parents per node. Hill Climbing accepts the best scoring neighbor until no improving move remains. Tabu Search adds a list of the $\tau$ most recent modifications that cannot be immediately reversed; $\tau =0$ recovers Hill Climbing exactly, and we use $\tau =15$. GES instead searches over Markov equivalence classes through a forward phase adding edges and a backward phase removing them, returning a completed partially directed graph converted to a single DAG for direct comparison, a genuinely different search paradigm and the reason we include it as the third algorithm.

3.4. Scoring Functions

Both scoring functions decompose additively across variables, $\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\left(G,D\right)={\displaystyle \sum _{i=1}^n}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\left(X_i,{\mathrm{P}\mathrm{a}}_G\left(X_i\right),D\right)$. BIC combines maximized log likelihood with a parameter count penalty,

$\begin{array}{r}\mathrm{B}\mathrm{I}\mathrm{C}\left(G,D\right)=\mathrm{l}\mathrm{o}\mathrm{g}P\left(D\mid G,\widehat{\theta }\right)-{\displaystyle \frac{d}{2}}\mathrm{l}\mathrm{o}\mathrm{g}N,\end{array}$with $d$ free parameters and $N$ observations. BDeu computes an exact marginal likelihood under a Dirichlet prior with equivalent sample size $\alpha$ (set to the conventional default of ten),

$\begin{array}{r}{\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}_{\mathrm{B}\mathrm{D}\mathrm{e}\mathrm{u}}\left(X_i,{\mathrm{P}\mathrm{a}}_G\left(X_i\right)\right)={\displaystyle \sum _{j=1}^{q_i}}\left[\mathrm{l}\mathrm{o}\mathrm{g}{\displaystyle \frac{\Gamma \left(\alpha /q_i\right)}{\Gamma \left(\alpha /q_i+N_{ij}\right)}}+{\displaystyle \sum _{k=1}^{r_i}}\mathrm{l}\mathrm{o}\mathrm{g}{\displaystyle \frac{\Gamma \left(\alpha /\left(q_i{r}_i\right)+N_{ijk}\right)}{\Gamma \left(\alpha /\left(q_i{r}_i\right)\right)}}\right],\end{array}$with $N_{ijk}$ the count of observations with $X_i$ in state $k$ and parents in configuration $j$. K2 uses the same Dirichlet-multinomial family as BDeu but fixes all pseudocounts at one regardless of parent configuration count, so its implied prior strength grows with the number of configurations, a structurally different, well established assumption.

3.5. Comparison, Bootstrap, and Sensitivity Checks

Structures are compared with the Structural Hamming Distance,

$\begin{array}{r}\mathrm{S}\mathrm{H}\mathrm{D}\left(G_1,G_2\right)=|{\tilde{E}}_1\backslash {\tilde{E}}_2|+|{\tilde{E}}_2\backslash {\tilde{E}}_1|+\left|\right\{e\in {\tilde{E}}_1\cap {\tilde{E}}_2:{\mathrm{d}\mathrm{i}\mathrm{r}}_1\left(e\right)\ne {\mathrm{d}\mathrm{i}\mathrm{r}}_2\left(e\right)\left\}\right|,\end{array}$over undirected skeletons $\tilde{E}$. To test whether algorithm agreement is stable under resampling, we generate $B=100$ bootstrap resamples per dataset and relearn all three algorithms on each, recording mean, median, and a 95% interval of the pairwise SHD for each of the three algorithm pairs, following the general logic of bootstrap edge reliability assessment used elsewhere (Prestat et al. 2013). We also pool all bootstrapped algorithm-pair SHDs against the full-sample scoring-function SHDs to summarize which source of disagreement is larger on average; because the scoring-function comparison is not itself bootstrapped in this study, this pooled comparison is reported with that caveat explicit (Section 6). We repeat the full pipeline at $k\in \{3,4,5\}$ on the reference dataset to test whether results depend on this choice.

4. Results

4.1. Datasets and Feature Selection

Table 1 summarizes the three datasets and the number of predictors each retained. Breast Cancer Wisconsin’s selected features are listed in Table 2; Pima Indians Diabetes retained all eight original predictors; Heart Disease retained nine of thirteen, dropping four (trestbps, chol, fbs, restecg) for redundancy with the retained set.

4.2. Search Algorithm Comparison Under BIC

Table 3 reports Hill Climbing, Tabu Search, and GES under BIC on all three datasets. Hill Climbing and Tabu Search are identical on every dataset (SHD = 0). GES matches both on Heart Disease exactly and differs by a single edge on Breast Cancer Wisconsin and Pima Diabetes, reaching the same BIC score in every case despite not directly following the same neighborhood-search trajectory.

4.3. Bootstrap Stability

Table 4 reports pairwise SHD across 100 bootstrap resamples per dataset for all three algorithm pairs. Hill Climbing and Tabu Search remained identical on essentially every resample (median SHD 0.0; mean 0.0 to 0.01). GES diverged from both considerably more once resampling is taken into account: median SHD 7.0 (Breast Cancer Wisconsin), 5.0 (Pima Diabetes), 3.0 (Heart Disease), compared with an SHD of 0 to 1 between the same pairs on the full sample (Table 3). Figure 1 shows per edge bootstrap frequencies for Hill Climbing and Tabu Search on Breast Cancer Wisconsin, identical for both.

4.4. Scoring Function Comparison

Table 5 compares BIC, BDeu, and K2 under Hill Climbing. BDeu matches BIC exactly on Pima Diabetes (SHD = 0) but differs by 7 and 8 edges on Breast Cancer Wisconsin and Heart Disease respectively; K2 differs from BIC on all three datasets (SHD 9 to 13), consistently producing denser networks. Figure 2 shows the BIC versus BDeu comparison for Breast Cancer Wisconsin. Figure 3 summarizes Structural Hamming Distance for both comparisons across all three datasets.

Table 6 pools all bootstrapped algorithm-pair SHDs (900 values) against all full-sample scoring-function SHDs (6 values). Algorithm differences have a lower median (3.0) than scoring function differences (8.5), though the latter are full-sample point estimates rather than bootstrapped (Section 6).

4.5. Sensitivity to the Redundancy Weight and the Discretization Level (Reference Dataset)

Table 7 reports the redundancy weight sweep on Breast Cancer Wisconsin: family diversity peaks between 1.3 and 1.5, where the selected feature set is identical. Table 8 reports the discretization level sweep: feature selection is stable (10 of 11 shared at levels 4 and 5 versus level 3), edge count decreases from 14 at level 3 to 11 at levels 4 and 5, and Hill Climbing and Tabu Search remain identical to each other (SHD = 0) at every tested level.

5. Discussion

Hill Climbing and Tabu Search agreed almost perfectly on all three datasets, both on the full sample and across 100 bootstrap resamples each (median SHD 0.0 throughout), consistent with a benchmark comparison where the two performed almost identically (Ribeiro and Maciel 2023). GES matched both closely on the full sample of every dataset (SHD 0 to 1) despite a structurally different search space, consistent with the equivalence class framing that motivates it (Chickering 1996); the bootstrap, however, shows this agreement does not extend to resampling stability, with median SHD between GES and either neighborhood based algorithm ranging from 3.0 to 7.0 across datasets. The three algorithms can therefore agree closely on one point estimate while differing considerably in sensitivity to sampling variation, a distinction a full-sample-only comparison would not reveal.

Scoring function disagreement was substantial on two of three datasets and absent on the third (Pima Diabetes, where BDeu matched BIC exactly). Pooling across all comparisons (Table 6), scoring function differences had a higher median SHD (8.5) than algorithm differences (3.0), consistent with reports that BIC recovers reference structures more accurately than BDeu (Liu et al. 2012; Liao et al. 2022) and that BDeu requires tuning of its equivalent sample size to remain competitive (Chua and Ong 2017). This pooled comparison mixes a bootstrapped distribution (900 values) with six full-sample point estimates, so it is suggestive rather than a matched statistical test. K2 diverged from BIC more than BDeu did on every dataset, consistent with K2’s fixed, configuration count independent pseudocount inflating effective prior strength as parent sets grow, favoring denser networks; a theoretical account of BDeu’s related tendency to add uninformative parents (Suzuki 2016) is directionally consistent with this pattern for both non-BIC scores.

Based on these three datasets, the scoring function was, on average, a larger source of structural disagreement than the search algorithm, but the algorithm comparison is not uniformly stable either: Hill Climbing and Tabu Search were interchangeable throughout, while GES’s agreement with them held on the full sample but not under resampling. Whether either pattern holds on further datasets is a question for replication.

6. Limitations and Future Work

This study used three datasets of moderate size, one discretization method, and three scoring functions; results should not be assumed to generalize without replication elsewhere, particularly at higher dimensionality (Prestat et al. 2013; Wang et al. 2023; Su et al. 2013). The redundancy weight and discretization level checks were run on the reference dataset only, for tractability. The scoring function comparison was not itself bootstrapped, so the pooled comparison against bootstrapped algorithm differences (Table 6) is suggestive rather than a matched statistical test. Quantile discretization discards within bin variation, and edge orientation under score based search does not imply causation. Future work could extend this comparison to further datasets, to genomic scale data, to additional scoring functions such as factorized normalized maximum likelihood (Liu et al. 2012), to stabilized search variants (Kitson and Constantinou 2025), and to a bootstrapped scoring function comparison placing both axes on equal statistical footing.

7. Conclusion

Across three healthcare datasets, Hill Climbing, Tabu Search, and GES recovered nearly identical structures under BIC on the full sample (SHD 0 to 1). A symmetric 100-resample bootstrap for all three algorithms showed this holds robustly for Hill Climbing versus Tabu Search (median SHD 0.0 on every dataset) but not for GES, which diverged more from both under resampling (median SHD 3.0 to 7.0) despite matching them on the full sample. Switching the scoring function changed the structure substantially on two of three datasets and not at all on the third; pooled across all comparisons, scoring function differences had a higher median SHD (8.5) than bootstrapped algorithm differences (3.0), though this mixes a bootstrapped distribution with full-sample point estimates. Two internal checks showed this pattern does not rest on an arbitrary choice of redundancy weight or discretization level. Within these three datasets, the scoring function was generally the larger source of disagreement, but algorithm-level agreement was itself uneven once resampling was considered; both findings require confirmation on further healthcare datasets before being treated as general.