Explainable AI for Stroke Prediction: A Hybrid Bayesian Network Approach

1. Introduction

Stroke is a leading cause of death and disability in Canada and worldwide [1]. It has been estimated that up to 80% of strokes are preventable through timely identification and targeted screening of risk factors, including hypertension, smoking, diabetes and other vascular risks [2]. In the United States, stroke mortality has declined by more than 70% over the past 50 years through advancements in prevention, acute care, and secondary risk reduction [2].

However, substantial disparities persist globally, which drives uneven access to screening, timely care, and health education [3]. For example, approximately 40% of patients presenting to emergency departments with acute stroke are unaware of warning signs and etiological risk factors [4].

The shift toward precision medicine has fundamentally reshaped the paradigms of disease risk stratification, moving away from “one-size-fits-all” screening toward data-driven, personalized risk assessment [5]. As established by Steyerberg et al. [6], the methodology of clinical prediction has expanded significantly beyond the foundational use of logistic regression and linear classifiers. While these models are the gold standard for transparency, the increasing availability of high-resolution electronic health records (EHR) has demanded more flexible algorithmic approaches. This evolution spans from rule-based expert systems and decision trees to probabilistic frameworks and sophisticated deep learning models [6,7,8]. While the increasing complexity of algorithmic models captures the intricate interplay among disease risk factors, it obscures the transparency of the underlying mathematical models. The schematic architecture of these models is shown in Supplementary material, Figure S1.

The simpler condition-based model, rule-based learning, relies on predefined sets of rules and represents a highly transparent approach. This transparency and logic of explicit conditions inspired the automated, hierarchical data splitting used in decision trees [9]. Similarly, K-nearest neighbors (KNN) and support vector machines (SVMs) solve nonlinearity by classifying patient risk mathematically via spatial distance and geometric hyperplanes [10,11]. Generalized additive models (GAMs) provide flexibility by isolating continuous variables (e.g., blood pressure) into standalone risk curves [12], whereas Bayesian networks use directed acyclic graphs to calculate explicit conditional probabilities between interdependent risk factors [6]. At the upper end of the model’s structural complexity, Ensembles (e.g., Random Forests and XGBoost) aggregate hundreds of individual models to generate a weighted consensus, and deep learning networks process high-dimensional variables through layers of hidden computational nodes [13].

Machine learning models, particularly neural networks and XGBoost, often achieve high predictive performance [14]. However, these models are frequently classified as “black-box” models, as the underlying reasoning of diagnostic decisions remains obscured from human interpretation [15,16]. This lack of transparency presents a critical barrier in high-stakes medical settings, where both healthcare providers and patients are often unable to discern the underlying biological rationale required to confidently trust and act upon the algorithmic output.

Conversely, traditional statistical approaches such as logistic and linear regression remain cornerstone tools in clinical risk prediction [17]. Within the formal taxonomy of explainable artificial intelligence (XAI), these models are classified as “transparent” [18], “white-box” [19] or “glass-box” [20] models, which characterize algorithms by their intrinsic interpretability.

Interpretability is formally defined as the degree to which a human can consistently predict the model’s result and comprehend the underlying cause of a decision [21]. It is the extent to which the internal mechanics of a model make sense to a human observer without requiring supplementary, post hoc explanations.

Arrieta et al. [18] characterize transparency through 3 hierarchical dimensions as listed below:

• Simulatability: A clinician can mentally or manually calculate a patient's stroke risk via a simple scoring system (e.g., adding points for age and hypertension).
• Decomposability: A physician can isolate and understand the specific weight (coefficient) assigned to a single variable, such as observing how a 10 mmHg increase in blood pressure explicitly shifts the risk output.
• Algorithmic transparency: The mathematical guarantee that a model will converge on a unique, optimal solution for a patient's risk profile, driven by a visible optimization process such as maximum likelihood estimation.

Expanding upon structural transparency, Doshi-Velez and Kim [22] propose a rigorous taxonomy for evaluating interpretability through three distinct task-driven categories.

• Application-grounded evaluations: This involves conducting experiments with domain experts within the context of real-world tasks, such as a neurologist evaluating the model’s output while diagnosing acute stroke patients in the emergency department.
• Human-grounded evaluations: Layperson or nonspecialist clinicians are tasked with identifying the most weighted risk factor from visual dashboards or charts.
• Functionally grounded evaluations: Researchers benchmarking a highly complex model (like a neural network) against mathematically transparent models (like a decision tree) without requiring live human experiments.

While these transparent architectures ensure clinical accountability, their intrinsic interpretability frequently comes at the expense of predictive accuracy [23]. Although these models may perform optimally within highly structured and well-defined datasets, they lack the mathematical flexibility required to map the widespread, nonlinear dependencies inherent in real-world clinical data [18]. Consequently, developing targeted diagnostic tools requires navigating a trade-off between maximizing performance and preserving structural transparency, as illustrated in Supplementary material, Figure S2.

The trade-off exhibits a Pareto frontier, where models inherently sacrifice mathematical visibility to achieve higher dimensionality. Arrieta et al. outlined two primary research trajectories for modern clinical diagnostics [18]. The first involves applying post hoc explainability techniques such as Shapley additive explanation (SHAP) [24] or local interpretable model-agnostic explanations (LIME) [25]. This is to estimate feature importance and approximate the logic of highly accurate, opaque models. Alternatively, researchers are developing hybrid modeling approaches that enhance the performance of interpretable models via advanced preprocessing techniques and hybrid models.

Lolak et al. [26] compared the performance and interpretability of logistic regression, Cox proportional hazard, Bayesian network, tree-augmented naïve Bayes (TAN), extreme gradient boosting (XGBoost), and explainable boosting machine (EBM) methods for predicting stroke risk. They included a retrospective cohort of 275,247 high-risk patients with discretized risk factors, including demographic, cardiovascular, and metabolic variables and medication use. Among these patients, 9659 (3.5%) experienced strokes. Multiple imputation by chained equation (MICE) was implemented for missing risk factors. The C statistics of XGBoost, EBM, TAN, logistic regression and the Bayesian network were 0.89, 0.87, 0.83, 0.80 and 0.80, respectively. However, XGBoost needs post hoc analysis to improve explainability. On the other hand, the explainability of these models has been less explored.

The optimal diagnostic model is dictated by the dataset’s architecture and the need for interpretability [15], as the “No Free Lunch” theorem establishes that no single algorithmic approach is universally superior across the entire dataset [27]. While deep learning excels at extracting features from massive, unstructured datasets, tree-based ensembles such as XGBoost frequently outperform neural networks on structured, tabular electronic health record (EHR) data [13,28]. Conversely, parsimonious algorithms such as logistic regression and Bayesian networks offer superior mathematical interpretability in cohorts characterized by limited sample sizes or significant data missingness [6]. Furthermore, selecting clinically relevant nonredundant features via effective feature selection not only mitigates the “curse of dimensionality” but also ensures that the resulting diagnostic model remains both mathematically sound and biologically plausible [6].

To bridge the gap between predictive performance and model interpretability, Bayesian networks (BNs) present a parsimonious alternative. On the accuracy‒interpretability continuum, BNs occupy a strategic intermediate position that bridges simpler transparent models and high-performing opaque architectures. The mathematical foundation of a Bayesian network (BN) is Baye’s theorem, which provides a rigorous framework for updating the probability of a clinical hypothesis ($H$, such as a stroke diagnosis) on the basis of observed evidence ($E$, such as specific patient risk factors):

$$P\left(H|E\right)={\displaystyle \frac{P\left(E|H\right)P\left(H\right)}{P\left(E\right)}}$$

The term $P\left(H|E\right)$ represents the posterior probability, which is the probability of the hypothesis given the evidence. In a clinical context, this allows for dynamic risk assessment where the prior probability $P\left(H\right)$, which is based on population prevalence, is redefined into a personalized risk score as individual patient data $\left(E\right)$ are integrated.

A Bayesian network extends this principle into a probabilistic graphical model that represents the joint probability distribution of a complex set of clinical variables. As illustrated in the sequential architecture of Figure 1, the model decomposes high-dimensional data into a product of local conditional distributions via the following theorem:

$$P\left(X_1,\dots .,X_n\right)={\displaystyle \prod _{i=1}^n}P\left(X_i\right|Parents\left(X_i\right))$$

The conditional probability table (CPT) serves as the primary parameterization of the network. Each CPT quantifies the probabilistic relationship between a "child" node and its "parent" nodes by explicitly listing the probability of every possible child state given all combinations of parent states. For example, a CPT for the node "Stroke" calculates specific risk probabilities on the basis of discrete parent states such as "Old Age" and "High Glucose". This modular structure ensures that the model remains decomposable, as individual tables can be isolated and audited by clinicians to verify that the algorithmic logic remains biologically plausible and aligned with clinical domain expertise.

Bayesian networks have been implemented in disease diagnosis and risk prediction, including for cardiovascular disease [29], atrial fibrillation relapse [30], and ALS diagnosis [31]. Suo et al. [32] developed a BN-based model using EHR data from 169,692 patients in Chinese hospitals to predict 3-year CHD risk. A weighted survival Bayesian network was developed by combining inverse probability censoring weights and a Cox model using demographic, laboratory and medical history variables as risk factors. The model achieved an AUC of 0.80 despite extensive missing data and imbalanced datasets.

Similarly, Park et al. [33] applied a BN to predict 3-month functional outcomes (mRS) and 1-year mortality in 2,432 ischemic stroke patients using 74 features (risk factors), including baseline NIHSS scores and age, attaining an AUC of 0.89 for 3-month functional outcomes and 0.89 for 1-year mortality. They also developed a Yonsei outcome inference system via their trained model.

Lolak et al. [26] demonstrated the utility of BNs in a retrospective cohort of 275,247 high-risk patients compared with other machine learning models. Although XGBoost achieved a superior C statistic of 0.89, the BN model demonstrated a competitive value of 0.80 while maintaining intrinsic interpretability. with numerous features.

Despite their utility in a large cohort and ability to capture many clinical features, BNs often lag in discriminatory performance. This predictive gap arises from inappropriate handling of missing data, inadequate preprocessing techniques and imbalanced targeted classes.

To address this limitation without sacrificing the intrinsic transparency of the BN model, our research aligns with the objective of increasing the predictive capability of BNs. By implementing rigorous data preprocessing techniques and utilizing a hybrid modeling structure to accommodate continuous variables, this study aims to close the accuracy gap.

We present a rigorous, eight-arm methodological ablation study evaluating the predictive architecture of Bayesian networks. We systematically isolate the effects of network topology (purely discretized vs. hybrid continuous nodes), imputation strategies (baseline median imputation vs. multiple imputation by chained equations [MICE]), and class balancing (unadjusted vs. synthetic minority oversampling technique [SMOTE]).

By isolating the effect of each methodological variation, the following hypothesis has been constructed:

Focus Area	Null Hypothesis (H0)	Alternate Hypothesis (H1)
1. Topology	Hybrid continuous topology does not significantly change predictive performance vs. pure discretization.	Hybrid continuous topology significantly changes predictive performance vs. pure discretization.
2. Imputation	MICE does not significantly improve performance metrics vs. baseline median imputation.	MICE significantly improves performance metrics vs. baseline median imputation.
3. Balancing	SMOTE does not significantly increase minority class sensitivity vs. unadjusted distributions.	SMOTE significantly increases minority class sensitivity vs. unadjusted distributions.
4. Synergy	The integration of continuous topology, MICE, and SMOTE does not significantly improve overall capacity vs. baseline.	The integration of continuous topology, MICE, and SMOTE significantly improves overall capacity vs. baseline.

2. Materials and Methods

Ethical approval and informed consent were not required for this study. We used an open-source Stroke prediction dataset, which is available at Kaggle [34]. The dataset has been widely adopted in machine learning research, with approximately 290,000 downloads, and is used in 103 research articles [35].

2.1. Overview of the Dataset

The dataset consists of 10 features along with the diagnosis of stroke. The dataset is highly imbalanced, with 4861 nonstroke patients and 249 stroke patients with a stroke incidence of 4.9%. The summary of the feature is shown in Supplementary material, Table S1.

The dataset features incorporate continuous, binary and categorical values. These features are broadly categorized in the following aspects:

• Sociodemographics: age, sex, residence type, marital status and occupation
• Medical history: Hypertension and heart disease
• Physiological measurements: BMI (body‒mass index) and average glucose level
• Lifestyle factor: Smoking status

Reshetnyak et al. demonstrated that cumulative sociodemographic factors serve as significant, independent predictors of incident stroke [36]. Similarly, the landmark INTERSTROKE study in the Lancet revealed that risk factors such as hypertension, cardiac history, BMI, glucose levels, and smoking status collectively account for approximately 90% of the population-attributable risk for stroke [37].

Continuous variables need to be converted to categorical variables for the Bayesian network. Age was categorized on the basis of the WHO’s standard age group (differences of 5 and 85+) [38]. Similarly, body mass index (BMI) has also been categorized on the basis of the WHO [39]: underweight (less than 18.4 kg/m2), normal weight (18.5 to 24.9 kg/m2), overweight (25 to 29.9 kg/m2) and obese (greater than 29.9 kg/m2). Diabetes was classified as present if the average glucose was ≥126 mg/dL, prediabetes if 100–125 mg/dL, and normal if <100 mg/dL [40]. This glucose-only approach ensures consistency in aggregated data while acknowledging limitations such as ignoring timing, symptoms, A1C, medications, and comorbidities that influence accurate diagnosis.

Alternately, we maintain these variables in their continuous form by employing a hybrid Bayesian network (HBN) structure [41]. To implement this hybrid approach, the following steps were taken:

• Age, BMI and average glucose were preserved as exact numeric values.
• The conditional distribution of continuous nodes was specified by a Gaussian function via conditional linear Gaussian distributions.
• We utilized structure and parameter learning algorithms explicitly designed for mixed-data networks to estimate outcomes accurately.

2.2. Study Design

This study employs a 2x2x2 factorial design to evaluate the impact of different data preparation techniques on model performance. Figure 1 illustrates the overall methodological framework, distinguishing between standard pipeline procedures and experimental variables. The baseline methodology (Model M1) establishes a fixed processing pipeline that is applied uniformly across all the models. To prevent data leakage, a standard training-validation test (60-20-20) split is executed prior to any data manipulation. Following the split, preprocessing choices are applied: median missing value imputation, feature discretization, and no class balancing.

2.3. Study Arms

To evaluate potential improvements over the baseline, three experimental variations were introduced at key methodological decision points (denoted by dashed lines in Figure 1):

• Imputation: Median vs. MICE
• Feature type: Discretized vs. Continuous
• Training Set Balancing: None vs. SMOTE

Based on the permutations of these three variables, the data were processed through eight distinct analytical study arms (Models M1 through M8), as detailed in Table 1. This configuration effectively isolates the individual and combined effects of MICE, continuous features, and SMOTE, allowing for a comprehensive comparison of how each preprocessing variation impacts overall predictive accuracy.

2.4. Structure Learning

The structure of the Bayesian networks (BNs) was learned via the hill-climbing (HC) heuristic search algorithm. To reduce computational overhead and ensure clinical validity, we implemented whitelisting and blacklisting techniques. The strongly associated nodes were whitelisted, and biologically impossible edges (such as blood pressure causing age) were blacklisted.

Different structural constraints were applied depending on whether the experimental arm utilized purely discrete or hybrid (mixed continuous-discrete) data. For hybrid arms, specialized whitelists and blacklists are required to accommodate the mathematical constraints of conditional Gaussian Bayesian networks (CGBNs), which restrict certain topological relationships (e.g., preventing continuous nodes from acting as parents to discrete nodes). To guarantee the robustness of the learned structures and avoid overfitting to the training set, model averaging by bootstrapping was employed. We generated 200 bootstrap iterations for each network, utilizing the resulting arc strengths to construct a robust consensus network.

2.5. Parameter Learning

Following structure learning, conditional probability tables (CPTs) and distributions were estimated from the training data. For purely discrete experimental arms (e.g., M1, M2), parameter learning was conducted via Bayesian estimation. To prevent zero-probability errors for rare feature combinations, an imaginary sample size (ISS) of 10 was applied, effectively acting as a smoothing prior. For hybrid networks containing continuous nodes (e.g., M3, M5), standard maximum likelihood estimation was used to fit the local probability distributions.

2.6. Probabilistic Inference

To generate stroke predictions on the evaluation sets, distinct inference strategies are required for the data types. For discrete networks, posterior probabilities were calculated rapidly via exact inference. Conversely, exact inference is computationally intractable for complex hybrid BNs. Therefore, for arms containing continuous variables, approximate inference was conducted via Monte Carlo likelihood weighting. For every individual patient in the evaluation set, 1,000 independent Monte Carlo simulations were run against the fitted hybrid network to derive a stable probability of the target event (stroke). Population-level counterfactual queries reported in the Results instead used logic sampling with 1,000,000 particles to obtain low-variance marginal estimates.

2.7. Model Evaluation and Statistical Analysis

Model performance was assessed by comparing the predicted probabilistic outputs against the true stroke outcomes in the blind hold-out test sets. Because stroke datasets are inherently imbalanced, a default classification threshold of 0.5 is often suboptimal. Therefore, for each model, a receiver operating characteristic (ROC) curve was generated, and the optimal classification threshold was determined dynamically via the Youden index (Youden's J statistic) [42]. This method selects the threshold that maximizes the sum of sensitivity and specificity, ensuring a balanced clinical trade-off between identifying true stroke cases and minimizing false alarms.

Using the optimal, model-specific thresholds, patients were classified into binary outcomes, and confusion matrices were constructed. The primary metrics extracted included the area under the ROC curve (AUC), sensitivity, specificity, and Youden's J statistics.

To ensure statistical rigor and evaluate the generalizability of the findings, 95% confidence intervals (CIs) for the AUC were calculated for every model. This was achieved via a robust bootstrap resampling procedure with 2,000 iterations per model. Probability calibration was quantified with the Brier score and decile reliability curves. The eight arms were ranked primarily by AUC and secondarily by Youden’s J. Pairwise differences were tested with DeLong’s test for correlated ROC curves (AUC) and McNemar’s test on the true-positive cohort (sensitivity). To confirm that the topology effect was not specific to a single partition, the entire pipeline was repeated over 30 independent stratified resamples, and the per-seed AUC differences between continuous and discrete arms were compared using paired t-tests and Wilcoxon signed-rank tests. Finally, logistic regression and an early-stopped gradient-boosted ensemble (XGBoost) were trained on identical splits as opaque benchmarks. Analyses were performed in R 4.5.2 (bnlearn 5.1) with a fixed global seed; package versions are recorded in the session information.

2.8. Software Version and System Specification

We implemented the “bnlearn” package version 5.1 in RStudio version 2025.09.2 Build 418 with R 4.5.2 to build our BN model. This study was conducted on a processor with an Intel Core Ultra 9 24 core and an NVIDIA GeForce RTX 5090 laptop GPU.

3. Results and Discussion

The distribution of each feature between each class is illustrated in Table 2. Similarly, the visualization of the feature distribution chart is shown in Table 2. The dataset is highly imbalanced, with only 4.9% of the population affected by stroke. The means and standard deviations are reported for continuous variables (age, average glucose level and BMI), whereas categorical variables are summarized using the frequency and percentage of the major category of each feature. P values were calculated via appropriate univariate statistical tests based on variable type, such as the chi-square test for categorical variables and the independent samples t-test for continuous variables.

From Table 2, we can observe several trends:

• Age: Stroke patients were older than nonstroke patients were (67.73 ± 12.73 vs 41.97 ± 22.29 years). Stroke cases are concentrated in the 50+ age brackets, peaking at approximately the 75--84 range.
• Hypertension: Hypertension was more common in the stroke group (26.5% vs 8.9%).
• Heart disease: Heart disease was more common among stroke patients (18.1% vs 4.7%).
• Marital status: The married population was more common in the stroke group (88.4% vs 64.5%).
• Work type: Most individuals in both groups work in the private sector (Private: 59.8% vs 5.1%). From the barchart, we can observe that self-employed individuals are relatively more common in the stroke group than in the baseline population.
• Residence type: Residence was similar between groups, with no meaningful difference (54.2% urban vs 50.6% urban).
• Glucose level: The average glucose level was greater in stroke patients (132.54 ± 61.92 vs 104.79 ± 43.85). According to the bar chart diagram, the proportion of patients in the stroke group was greater in the prediabetes and diabetes categories than in the nonstroke baseline group.
• BMI was greater in the stroke group (30.47 ± 6.33 vs 28.82 ± 7.91). Categorically, stroke patients are clustered in the “overweight” and “obese” brackets, with few falling into the “normal weight” category compared with nonstroke patients. Additionally, 201 patients were missing BMI data.
• Smoking status: While “never smoked” was the most common known status for both groups, the proportion of “formerly smoked” patients was greater in the stroke group. Additionally, a significant portion of the dataset (1544 patients overall) has an “unknown” smoking status.
• Overall, the stroke group presented a greater burden of vascular and metabolic risk factors.

3.1. Significance Table

An 11x11 significance matrix table is constructed, as shown in Supplementary material,Table S4. Owing to mixed data types (continuous, binary, and nominal categorical), we implemented the following statistical tests tailored to each pair of variables:

• Continuous vs continuous: We used Spearman’s rank correlation, which is a nonparametric test.
• Categorical vs Categorical: A chi-square test of independence was applied.
• Continuous vs Categorical: We implemented the Kruskal‒Wallis test, which tests the median of continuous variables across the groups of categorical variables.

The effect size (Spearman’s $\rho$ for continuous variables, Cramer’s V for categorical variables, and $\mathsf{\eta }$(Eta) derived from the Kruskal‒Wallis test for mixed pairs) and the significance are illustrated.

The resulting matrix reveals distinct and clinically significant clusters. Age was the most notable driver of both stroke (0.25***) and other lifestyle and clinical variables: marital status (0.66***), work type (0.64***), BMI (0.38***) and smoking status (0.38***). Similarly, hypertension (0.13***) and heart disease (0.13***) emerged as the next strongest predictors of stroke. Furthermore, both residence type and sex had near-zero effect sizes on stroke (0.01 each) and other variables. Consequently, age and cardiovascular indicators (hypertension and heart disease) related to stroke nodes were whitelisted for predictive modeling, whereas noncontributing factors such as residence type and sex were blacklisted.

3.2. Directed Acyclic Graph (DAG)

Based on performance metrics on Table 5 and Figure 2, Model 6 with MICE imputation with continuous variables is the best performing model. Further, bootstrap arc-strength and thickness is shown in Supplementary material, Table S4 and Figure S4.

Stroke is directly influenced by three discrete parent nodes, namely, heart_disease, ever_married and hypertension, along with two continuous nodes, namely, age and avg_glucose_level, as shown in Figure 2. Visually, the directed acyclic graph (DAG) depicts the discrete stroke node as a parent to continuous nodes such as age and avg_glucose_level. While physiologically counterintuitive, this is a known mathematical constraint of conditional Gaussian Bayesian networks (CGBNs). In the standard CGBN architecture, continuous variables cannot act as parent nodes to discrete variables. Therefore, to capture the strong statistical correlations between age, glucose, and stroke without violating the network's strict probabilistic rules, the structure learning algorithm inverts the arc direction. This ensures that the dependency is mathematically preserved and factored into the joint probability distribution, even if the graphical representation reverses the chronological causality. We can also observe that the residence_type node is completely isolated with no significant association with any of the nodes, whereas sex contributes to heart_disease.

3.3. Conditional Probability Table (CPT)

The conditional probabilities of all the nodes associated directly with the stroke are illustrated in Table 3 and Table 4.

The generative nature of our hybrid Bayesian network offers transparency, supporting the “white-box” paradigm. The exact probabilistic weights assigned to patient features are illustrated in Table 3. The Conditional Probability Table (CPT) for stroke illustrates how compounding risks dynamically escalate target probabilities. For example, a baseline patient (no heart disease, no hypertension, never married) has a minimal stroke probability of 1.19%. Introducing heart disease alone elevates this risk to 15.38%. Notably, the model captures synergistic multivariate interactions, demonstrating a severe probability spike to 75.00% when both heart disease and hypertension present simultaneously in unmarried individuals. Notably, this cell is supported by only four training patients (three of whom experienced a stroke), which accounts for the implausibly large swing relative to neighbouring cells. While the introduction of marriage decreases the localized probability for this specific comorbidity cluster to 13.33%, such extreme fluctuations within hyperspecific patient intersections likely reflect data sparsity rather than a definitive physiological rule, highlighting the need for cautious interpretation.

Table 4 highlights the mechanical advantage of retaining continuous architectures when the discrete topology is used. Rather than discretizing data into rigid bins, the network dynamically calculates expected values (e.g., average glucose level) via localized regression coefficients conditioned on discrete parent states. The expected continuous distribution is defined by the following formula:

$$ExpectedValue={\beta }_0+\left({\beta }_{1}xAge\right)+\left({\beta }_{2}xBMI\right)$$

These findings demonstrate that physiological relationships shift dynamically under varying clinical conditions. For example, the predictive weight of BMI (${\beta }_2$) is negligible at 0.016 for unmarried, stroke-free patients but multiplies substantially to 2.500 for married patients with a history of stroke. Concurrently, the standard deviation (σ) expands from 31.90 to 59.87, accurately reflecting the increased metabolic variance typically observed in poststroke populations. This localized parameterization confirms that hybrid networks successfully preserve natural physiological variance, avoiding the structural information loss inherent in arbitrary categorization.

3.4. Inference

Beyond static parameter evaluation, the primary clinical advantage of our generative framework lies in its capacity for dynamic probabilistic inference and counterfactual simulation. To illustrate this deployability, we conducted forward inference queries via logic sampling (n = 1,000,000 simulations) on the optimal hybrid network. For a baseline clinical persona (age ≤ 50, normotensive, and free of heart disease), the network yielded a nominal localized stroke probability of only 0.74%. Conversely, querying a high-risk clinical profile (age ≥ 65, hypertensive, positive for heart disease, and average glucose levels ≥ 150 mg/dL) significantly elevated the predicted probability to 43.18%.

Crucially, the generative architecture allows clinicians to simulate targeted medical interventions dynamically. If this high-risk patient's glucose levels were therapeutically managed to a normoglycemic state (≤100 mg/dL)—holding all other demographic and physiological variables constant—the network instantly recalibrated the stroke probability down to 31.55%. This capacity to transition from static black-box classification to interactive risk simulation underscores the profound clinical utility of the hybrid Bayesian architecture. It functions as a mathematically transparent decision-support tool capable of not only stratifying baseline risk but also quantifying the expected absolute risk reduction of specific preventative treatments. These queries are presented to demonstrate the inferential capability of the framework; because they derive from a single open-access dataset of uncertain provenance, the specific magnitudes should be regarded as illustrative of the method rather than as validated clinical effect estimates.

3.5. Performance Metrics

The performance metrics for each experimental arm with the baseline are shown in Table 5. Similarly, the ROC curve is illustrated in Figure 3. All the experimental groups were sorted based on the AUC followed by Youden index (J).

Although M5 (MICE + continuous + SMOTE) attained the highest Youden's J (0.488), the arms are ordered by AUC because it is threshold-independent; M5's apparent sensitivity advantage reflects a low operating threshold rather than superior discrimination, and it carries the worst calibration of any arm (Brier 0.182).

3.6. Hypothesis Testing

We conducted DeLong’s test [43] and McNemar’s test [44] to evaluate our research hypotheses. Delong’s test compares the overall predictive power (area under the curve, AUC) of correlated ROC curves evaluated on the same test cohort. McNemar’s tests were applied to the true-positive cohort to determine significant differences in minority class sensitivity.

Hypothesis 1. 

Topology (continuous vs. discrete)

To evaluate Hypothesis 1, DeLong’s and McNemar’s tests assessed overall discrimination and minority-class sensitivity, respectively. In the primary comparison against the baseline, continuous topology (M3) significantly improved overall discrimination (AUC 0.803 vs 0.718; p = 0.022) and true-positive stroke detection (McNemar sensitivity, p = 0.010) relative to the discrete baseline (M1). In the advanced-imputation framework, the continuous MICE model (M6) again significantly improved AUC over the discrete MICE model (M2) (p = 0.011), with a concordant but non-significant trend in sensitivity (p = 0.099). The null hypothesis is therefore rejected for discrimination in both imputation settings, establishing that continuous topology is superior to discretized features for clinical stroke prediction.

Hypothesis 2. 

Imputation (MICE vs.Median)

With respect to Hypothesis 2, the performance of the model with MICE applied to discrete variables produced no significant change in AUC (M2 0.720 vs M1 0.718). DeLong’s test revealed no statistically significant difference in predictive performance when MICE imputation was introduced to either the discrete topology (M2 vs M1, p = 0.708) or the continuous topology (M6 vs M3, p = 0.189). Therefore, we fail to reject the null hypothesis.

Hypothesis 3. 

Balancing (SMOTE vs.Unadjusted)

For Hypothesis 3, we fail to reject the null hypothesis. McNemar’s test on the true-positive cohort confirmed that SMOTE did not significantly increase minority-class sensitivity in any architectural variant: against the unadjusted baseline (M4 vs M1, p = 0.182), within the discrete MICE framework (M8 vs M2, p = 0.228), or within the continuous MICE framework (M5 vs M6, p = 0.450). Furthermore, SMOTE did not improve discrimination (M4 vs M1 AUC, p = 0.060) and roughly tripled the Brier score (from ~0.05 to ~0.18), indicating that synthetic oversampling degrades probability calibration without any compensating gain in sensitivity.

Hypothesis 4. 

Synergy (combinations vs. baseline)

Hypothesis 4 assessed the synergistic impact of combining methodologies against the strict unadjusted baseline (M1). Integrating all techniques (M5: continuous + MICE + SMOTE) failed to achieve a significant improvement over the baseline (Z = 1.77, p = 0.077), again reflecting the detrimental effect of synthetic balancing. In contrast, once SMOTE was ablated, the optimal configuration of continuous topology and MICE imputation (M6) achieved a highly significant improvement over the baseline (Z = 2.51, p = 0.012). Thus, the targeted synergy of continuous structuring—not the indiscriminate stacking of techniques—drives model enhancement.

As seen in Figure 4, the unbalanced models (M1, M3, M6) track the diagonal at low predicted risk, whereas the SMOTE model (M5) departs markedly, assigning predicted probabilities of 75–90% to patients whose observed stroke rate is ≈16%.

3.7. Robustness Across Resampling

To confirm that the topology advantage was not an artefact of a single partition, we repeated the entire pipeline across 30 independent stratified resamples. Continuous arms were consistently more discriminative than their discrete counterparts (M3 0.821 ± 0.028 vs M1 0.770 ± 0.039; M6 0.820 ± 0.028 vs M2 0.769 ± 0.040). The paired difference was highly significant for both the median-imputed pair (mean ΔAUC = 0.052; paired t = 8.14, p = 7.3 × 10⁻⁹; Wilcoxon p = 5.2 × 10⁻⁸) and the MICE-imputed pair (mean ΔAUC = 0.052; paired t = 8.19, p = 6.4 × 10⁻⁹). The continuous models also exhibited lower variance across resamples, indicating that preserving physiological variance yields discrimination that is not only higher but more stable (as shown in Figure 4; per-seed AUCs in Supplementary material, Table S3).

3.8. Calibration

Because a probabilistic decision-support tool must report reliable probabilities, we assessed calibration with the Brier score and reliability curves. The four unbalanced models were well calibrated (Brier 0.051–0.056), whereas every SMOTE-balanced model deteriorated roughly threefold (Brier 0.173–0.184). The reliability curves expose the mechanism: SMOTE models assigned predicted risks of 75–90% to patients whose observed stroke frequency was only about 16%, a systematic over-estimation, while the unbalanced models tracked the identity line. Synthetic oversampling therefore does not merely fail to improve the model; it corrupts the posterior probabilities that give a Bayesian network its clinical value (as shown in Table 6; full per-model Brier scores in Supplementary material, Table S2).

3.9. Comparison with Black-Box Models

To position the interpretable network against black-box alternative, an early-stopped gradient-boosted ensemble (XGBoost) were trained on identical splits and the continuous feature set. Similarly, we also trained logistic regression on the identical splits and continuous feature set. Discrimination was comparable across all three (Bayesian network 0.812; logistic regression 0.819; XGBoost 0.822). The interpretable network thus reached discrimination within approximately one AUC point of a tuned black-box ensemble while retaining full probabilistic transparency and counterfactual inference, directly addressing the accuracy–interpretability trade-off that motivates this study.

Findings from our study were compared with those in the literature: (paragraphs)

• We observed a statistically significant improvement in predictive power and stroke detection using continuous features compared with discretization. This finding supports the foundational statistical principle highlighted by Steyerberg [6], who cautions the categorization of continuous predictors. It also aligns with Murphy’s mathematical proofs regarding information retention in hybrid Bayesian networks (BNs) [41]. However, continuous models might sacrifice “risk bins” (e.g., diabetic vs. prediabetic patients), which clinicians naturally use in high-stake environments.
• We observed no statistically significant advantage in predictive performance when multivariate imputation (MICE) was used over simple median substitution. While MICE is widely recommended and utilized in the broader machine learning literature to handle missing healthcare data [7,8], its necessity in native probabilistic models appears less absolute. Interestingly, our highest-performing optimal model (M6) utilized MICE alongside continuous variables, suggesting that advanced imputation might offer subtle synergistic benefits within a continuous framework. Furthermore, these results must be interpreted with caution; in our specific dataset, missing values were primarily restricted to a single variable (BMI) and were likely missing completely at random (MCAR). Future studies should evaluate Bayesian missing data handling across datasets with heavier "missing not at random" (MNAR) patterns to establish exactly when computational imputation becomes strictly necessary over simple marginalization.
• Techniques such as the synthetic minority oversampling technique (SMOTE) have been widely adopted across machine learning methods to handle class imbalance and force algorithms to recognize minority classes [13,28]. However, our study revealed that the application of SMOTE failed to increase minority class sensitivity and actively degraded the overall performance of our models. This finding strongly resonates with emerging critiques (such as van den Goorbergh), which argue that synthetic oversampling distorts calibration and true prior probabilities. In a generative model such as a Bayesian network, synthesizing artificial patients forces "impossible" mathematical scenarios that might corrupt the native conditional probability tables. Consequently, future research addressing imbalanced clinical datasets might explore alternative balancing frameworks specifically designed for probabilistic architectures. For example, approaches such as the Bayesian network-based oversampling method (BOSME) [45] integrate Bayesian principles directly into the oversampling process for indirect cost-sensitive learning.
• We observed that integrating all advanced methodologies simultaneously (Continuous + MICE + SMOTE) yielded marginally significant results, whereas targeted synergy (combining continuous structuring with baseline robustness while dropping SMOTE) yielded highly significant improvements over the baseline. This observation aligns closely with the explainable AI (XAI) literature championed by Rudin [15], who cautions against overly complex algorithmic pipelines in high-stakes medical decisions. Future clinical modeling efforts may benefit from a more targeted approach, focusing on the careful validation of architectural choices that mathematically align with the natural distribution of the clinical data.

3.10. Limitations

Although the dataset widely used in machine learning models lacks verified provenances, its potential synthetic nature may limit the generalizability of our findings to real-world EHR data [35]. Future work should validate the model on prospectively collected data to confirm its performance in the hospital setting. Accordingly, the counterfactual and absolute-risk values reported here illustrate the model’s inferential mechanics and should not be read as clinical recommendations; external and prospective validation on verified electronic health records is required before any clinical deployment.

A known limitation of Conditional Gaussian Bayesian Networks (CGBNs) is the mathematical constraint prohibiting continuous variables from acting as parents to discrete nodes. Consequently, arcs between continuous predictors (Age, Glucose) and Stroke were directed outward from the target variable. While this inverts biological temporality, it mathematically preserves the joint probability distribution required for robust diagnostic and interventional inference.

4. Conclusions

This study demonstrated that preserving continuous physiological data within a hybrid conditional Gaussian Bayesian network significantly outperforms traditional data discretization for clinical stroke prediction. Our optimal synergistic architecture achieved a statistically significant improvement in predictive discrimination (AUC 0.812 vs. 0.718) over the unadjusted discrete baseline, challenging the common practice of categorizing metabolic variables and proving that preserving natural variance is critical for accurate risk stratification. Crucially, our findings caution against the indiscriminate application of standard machine learning preprocessing in generative architectures. Synthetic minority oversampling (SMOTE) failed to improve sensitivity and actively degraded overall performance, indicating that forcing artificial data points into mathematical models distorts native probability distributions.

The primary clinical value of this optimized hybrid topology lies in its utility as a transparent, "white-box" decision support system capable of dynamic probabilistic inference. These models excel at identifying highly specific at-risk populations; just as in our case where the network detected a severe 75.00% probability spike in a specific comorbid subpopulation, clinicians can pinpoint and quantify localized vulnerabilities. Furthermore, the model empowers users to simulate counterfactual therapeutic interventions, such as quantifying an absolute risk reduction from 43.18% to 31.55% following targeted glycemic management. However, a primary limitation of this study is its reliance on open-source data from unverified provenances. Pending future validation with prospectively collected, real-world electronic health records, this study establishes that methodologically conservative, continuous Bayesian networks offer a highly interpretable and mathematically sound framework for real-time diagnostic support.