Hybrid Bayesian Network for Interpretable Stroke Prediction — Final Analysis Pipeline

00 Setup & reproducibility

library(tidyverse)
library(bnlearn)
library(caret)
library(pROC)
library(mice)
library(recipes)
library(themis)

# One global seed governs the whole single-split pipeline.
GLOBAL_SEED <- 42
set.seed(GLOBAL_SEED)

01 Load & clean data

stroke_data <- read.csv("healthcare-dataset-stroke-data.csv")

stroke_data <- stroke_data %>%
  filter(gender != "Other") %>%      # single ambiguous record
  droplevels() %>%
  mutate(
    bmi            = as.numeric(as.character(bmi)),   # "N/A" strings -> NA
    stroke         = as.factor(stroke),
    hypertension   = as.factor(hypertension),
    heart_disease  = as.factor(heart_disease),
    gender         = as.factor(gender),
    ever_married   = as.factor(ever_married),
    work_type      = as.factor(work_type),
    Residence_type = as.factor(Residence_type),
    smoking_status = as.factor(smoking_status)
  )

cat(sprintf("N = %d patients | stroke prevalence = %.1f%%\n",
            nrow(stroke_data), 100 * mean(stroke_data$stroke == "1")))

N = 5109 patients | stroke prevalence = 4.9%

02 Exploratory data analysis (optional, for figures)

# Set eval=TRUE to regenerate the EDA figures used in the manuscript.
# (Bar charts for categorical variables, box plots for continuous variables.)
# --- kept from the original notebook; omitted here for brevity of the run ---

03 Stratified 60/20/20 split

# Stratified split that preserves the (rare) stroke prevalence in each fold.
# A `seed` argument is exposed so the SAME function can drive the repeated-split
# robustness analysis in Segment 14.
split_stroke <- function(data, seed = GLOBAL_SEED, p_train = 0.6, p_val = 0.2) {
  set.seed(seed)
  s0 <- data %>% filter(stroke == "0"); s1 <- data %>% filter(stroke == "1")
  n0 <- nrow(s0); n1 <- nrow(s1)

  tr0 <- sample(n0, round(n0 * p_train)); tr1 <- sample(n1, round(n1 * p_train))
  va0 <- sample(setdiff(seq_len(n0), tr0), round(n0 * p_val))
  va1 <- sample(setdiff(seq_len(n1), tr1), round(n1 * p_val))

  list(
    train = bind_rows(s0[tr0, ], s1[tr1, ]),
    val   = bind_rows(s0[va0, ], s1[va1, ]),
    test  = bind_rows(s0[setdiff(seq_len(n0), c(tr0, va0)), ],
                      s1[setdiff(seq_len(n1), c(tr1, va1)), ])
  )
}

splits   <- split_stroke(stroke_data, seed = GLOBAL_SEED)
trainset <- splits$train; valset <- splits$val; testset <- splits$test

for (nm in c("trainset", "valset", "testset")) {
  d <- get(nm); p <- prop.table(table(d$stroke))
  cat(sprintf("%-9s n=%4d | stroke=%.3f\n", nm, nrow(d), p["1"]))
}

trainset  n=3065 | stroke=0.049
valset    n=1022 | stroke=0.049
testset   n=1022 | stroke=0.049

04 Preprocessing helpers

# --- (a) Discretising imputer: median BMI from TRAIN only, then WHO/CDC/ADA bins
impute_baseline <- function(data, train_median_bmi) {
  data %>%
    mutate(bmi = as.numeric(as.character(bmi)),
           bmi = ifelse(is.na(bmi), train_median_bmi, bmi)) %>%
    mutate(
      age = cut(age, breaks = c(seq(0, 85, by = 5), Inf), right = FALSE,
                labels = c("0-4","5-9","10-14","15-19","20-24","25-29","30-34",
                           "35-39","40-44","45-49","50-54","55-59","60-64","65-69",
                           "70-74","75-79","80-84","85+")),
      bmi = cut(bmi, breaks = c(0, 18.5, 25, 30, Inf), right = FALSE,
                labels = c("Underweight","Normal","Overweight","Obese")),
      avg_glucose_level = cut(avg_glucose_level, breaks = c(0, 100, 126, Inf),
                right = FALSE, labels = c("Normal","Prediabetes","Diabetes"))
    ) %>%
    mutate(across(c(gender,hypertension,heart_disease,ever_married,
                    work_type,Residence_type,smoking_status,stroke), as.factor)) %>%
    droplevels() %>% select(-any_of("id"))
}

# --- (b) Continuous/hybrid prep: keep age/bmi/glucose numeric
prep_continuous <- function(df, train_median_bmi) {
  df %>%
    mutate(bmi = as.numeric(as.character(bmi)),
           bmi = ifelse(is.na(bmi), train_median_bmi, bmi),
           avg_glucose_level = as.numeric(avg_glucose_level),
           age = as.numeric(age)) %>%
    mutate(across(c(gender,hypertension,heart_disease,ever_married,
                    work_type,Residence_type,smoking_status,stroke), as.factor)) %>%
    select(-any_of("id"))
}

# --- (c) SMOTENC synthetic oversampling for mixed (categorical+continuous) data
apply_smotenc <- function(df, target_var = "stroke") {
  df <- droplevels(as.data.frame(df))
  rec <- recipe(formula(paste(target_var, "~ .")), data = df) %>%
    step_smotenc(all_outcomes(), over_ratio = 1, seed = GLOBAL_SEED) %>% prep()
  out <- juice(rec) %>%
    mutate(across(where(is.numeric), as.numeric),
           across(where(is.integer), as.numeric),
           across(where(is.character), as.factor))
  droplevels(as.data.frame(out))
}

05 Build the 8 ablation arms (M1–M8)

# Leakage-safe MICE: combine the folds but `ignore` everything except TRAIN rows,
# so imputation rules are learned from training data only.
trainset$split <- "train"; valset$split <- "val"; testset$split <- "test"
combined <- bind_rows(trainset, valset, testset) %>% select(-any_of("id"))
mice_safe <- suppressWarnings(
  mice(combined, m = 1, ignore = combined$split != "train",
       maxit = 5, seed = GLOBAL_SEED, printFlag = FALSE))
done <- complete(mice_safe, 1)
train_m2_raw <- done %>% filter(split=="train") %>% select(-split)
val_m2_raw   <- done %>% filter(split=="val")   %>% select(-split)
test_m2_raw  <- done %>% filter(split=="test")  %>% select(-split)
trainset <- trainset %>% select(-split)
valset   <- valset   %>% select(-split)
testset  <- testset  %>% select(-split)

med <- median(as.numeric(as.character(trainset$bmi)), na.rm = TRUE)  # TRAIN median

# Factorial 2x2x2: {discrete|continuous} x {median|MICE} x {none|SMOTE}
train_m1 <- impute_baseline(trainset, med); val_m1 <- impute_baseline(valset, med); test_m1 <- impute_baseline(testset, med)   # discrete/median/none
train_m2 <- impute_baseline(train_m2_raw, med); val_m2 <- impute_baseline(val_m2_raw, med); test_m2 <- impute_baseline(test_m2_raw, med) # discrete/MICE/none
train_m3 <- prep_continuous(trainset, med); val_m3 <- prep_continuous(valset, med); test_m3 <- prep_continuous(testset, med)   # continuous/median/none

train_m4 <- apply_smotenc(train_m1); val_m4 <- val_m1; test_m4 <- test_m1                                                     # discrete/median/SMOTE
train_m5_base <- prep_continuous(train_m2_raw, med)
train_m5 <- apply_smotenc(train_m5_base); val_m5 <- prep_continuous(val_m2_raw, med); test_m5 <- prep_continuous(test_m2_raw, med) # continuous/MICE/SMOTE
train_m6 <- train_m5_base; val_m6 <- val_m5; test_m6 <- test_m5                                                               # continuous/MICE/none
train_m7 <- apply_smotenc(train_m3); val_m7 <- val_m3; test_m7 <- test_m3                                                     # continuous/median/SMOTE
train_m8 <- apply_smotenc(train_m2); val_m8 <- val_m1; test_m8 <- test_m1                                                     # discrete/MICE/SMOTE
cat("Arms M1-M8 built.\n")

Arms M1-M8 built.

06 Feature association screen

# Cramer's V (categorical), point-biserial (continuous) vs stroke; used to
# justify which edges are whitelisted/blacklisted in Segment 07.
suppressPackageStartupMessages({library(vcd)})
cat_vars <- c("gender","hypertension","heart_disease","ever_married",
              "work_type","Residence_type","smoking_status")
num_vars <- c("age","avg_glucose_level","bmi")
cat_tbl <- map_dfr(cat_vars, function(v){
  tab <- table(stroke_data[[v]], stroke_data$stroke)
  data.frame(Variable=v, Assoc=round(assocstats(tab)$cramer,4),
             P=chisq.test(tab)$p.value)})
num_tbl <- map_dfr(num_vars, function(v){
  data.frame(Variable=v,
             Assoc=round(cor(stroke_data[[v]], as.numeric(as.character(stroke_data$stroke)),
                             use="complete.obs"),4),
             P=t.test(stroke_data[[v]]~stroke_data$stroke)$p.value)})
print(bind_rows(cat_tbl, num_tbl) %>% arrange(desc(abs(Assoc))))

            Variable  Assoc            P
1                age 0.2452 2.175773e-95
2      heart_disease 0.1349 2.120831e-21
3  avg_glucose_level 0.1320 2.373124e-11
4       hypertension 0.1279 1.688936e-19
5       ever_married 0.1083 1.686286e-14
6          work_type 0.0981 5.409035e-10
7     smoking_status 0.0756 2.007704e-06
8                bmi 0.0423 3.377378e-04
9     Residence_type 0.0154 2.998252e-01
10            gender 0.0091 5.598278e-01

07 Structural constraints (whitelist / blacklist)

all_nodes <- c("age","gender","hypertension","heart_disease","ever_married",
               "work_type","Residence_type","avg_glucose_level","bmi",
               "smoking_status","stroke")
discrete_nodes <- setdiff(all_nodes, c("age","avg_glucose_level","bmi","stroke"))

# Discrete models: clinically supported predictors point INTO stroke.
wl_discrete <- matrix(c("age","stroke","heart_disease","stroke",
                        "hypertension","stroke","avg_glucose_level","stroke"),
                      ncol=2, byrow=TRUE, dimnames=list(NULL,c("from","to")))
bl_discrete <- bind_rows(
  expand.grid(from=all_nodes, to=c("age","gender"), stringsAsFactors=FALSE) %>% filter(from!=to),
  expand.grid(from="stroke",  to=all_nodes, stringsAsFactors=FALSE) %>% filter(to!="stroke"),
  data.frame(from=c("gender","Residence_type"), to=c("stroke","stroke"))
) %>% as.matrix(); colnames(bl_discrete) <- c("from","to")

# Hybrid models: CGBN forbids continuous->discrete, so age/glucose arcs reverse.
wl_hybrid <- matrix(c("stroke","age","heart_disease","stroke",
                      "hypertension","stroke","stroke","avg_glucose_level"),
                    ncol=2, byrow=TRUE, dimnames=list(NULL,c("from","to")))
bl_hybrid <- bind_rows(
  expand.grid(from=all_nodes, to="gender", stringsAsFactors=FALSE) %>% filter(from!=to),
  expand.grid(from="stroke",  to=discrete_nodes, stringsAsFactors=FALSE) %>% filter(to!="stroke"),
  data.frame(from=c("gender","Residence_type"), to=c("stroke","stroke"))
) %>% as.matrix(); colnames(bl_hybrid) <- c("from","to")
cat("Whitelists/blacklists ready.\n")

Whitelists/blacklists ready.

08 Core evaluation function

# Trains one BN arm, tunes the threshold on VALIDATION, scores a blind TEST set.
# Structure learning: hill-climbing (HC). Set boot_R>0 for arc-strength bootstrap.
evaluate_variant <- function(train_df, val_df, test_df, model_name, boot_R = 200) {
  is_disc <- all(sapply(train_df, is.factor))
  wl <- if (is_disc) wl_discrete else wl_hybrid
  bl <- if (is_disc) bl_discrete else bl_hybrid

  dag <- hc(train_df, whitelist = wl, blacklist = bl)
  arc_str <- if (boot_R > 0)
    boot.strength(train_df, R = boot_R, algorithm = "hc",
                  algorithm.args = list(whitelist = wl, blacklist = bl)) else NULL
  fitted <- if (is_disc) bn.fit(dag, train_df, method = "bayes", iss = 10) else bn.fit(dag, train_df)

  # Posterior P(stroke=1): exact for discrete, likelihood-weighting for hybrid.
  get_probs <- function(td) {
    if (is_disc) {
      pp <- predict(fitted, node="stroke", data=td, method="bayes-lw", prob=TRUE)
      return(attr(pp,"prob")["1",])
    }
    ev_cols <- setdiff(names(td), "stroke")
    vapply(seq_len(nrow(td)), function(i){
      p <- suppressWarnings(cpquery(fitted, event=(stroke=="1"),
             evidence=as.list(td[i, ev_cols]), method="lw", n=1000))
      if (is.na(p)) 0 else p }, numeric(1))
  }

  val_probs  <- get_probs(val_df)
  roc_val    <- roc(val_df$stroke, val_probs, levels=c("0","1"), quiet=TRUE)
  opt_thresh <- coords(roc_val, "best", ret="threshold", best.method="youden")$threshold[1]
  if (is.na(opt_thresh)) opt_thresh <- 0.5

  test_probs <- get_probs(test_df)
  roc_test   <- roc(test_df$stroke, test_probs, levels=c("0","1"), quiet=TRUE)
  test_preds <- factor(ifelse(test_probs >= opt_thresh, "1", "0"), levels=c("0","1"))
  cm <- confusionMatrix(test_preds, test_df$stroke, positive="1")

  list(
    results = data.frame(Model=model_name, Val_Threshold=round(opt_thresh,4),
      Test_AUC=round(as.numeric(auc(roc_test)),4),
      Test_Sensitivity=round(cm$byClass["Sensitivity"],4),
      Test_Specificity=round(cm$byClass["Specificity"],4),
      Test_Youden_J=round(cm$byClass["Sensitivity"]+cm$byClass["Specificity"]-1,4),
      Test_F1=round(as.numeric(cm$byClass["F1"]),4)),
    roc=roc_test, fit=fitted, structure=dag, arc_strength=arc_str,
    confusion_matrix=cm, test_probs=test_probs)
}

09 MAIN RESULT — train & rank the 8 arms

experiments <- list(
  list(train=train_m1,val=val_m1,test=test_m1,name="M1: Baseline"),
  list(train=train_m2,val=val_m2,test=test_m2,name="M2: MICE + Discretize"),
  list(train=train_m3,val=val_m3,test=test_m3,name="M3: Continuous Nodes"),
  list(train=train_m4,val=val_m4,test=test_m4,name="M4: SMOTE on Baseline"),
  list(train=train_m5,val=val_m5,test=test_m5,name="M5: ALL (MICE + Cont + SMOTE)"),
  list(train=train_m6,val=val_m6,test=test_m6,name="M6: MICE + Continuous"),
  list(train=train_m7,val=val_m7,test=test_m7,name="M7: Continuous + SMOTE"),
  list(train=train_m8,val=val_m8,test=test_m8,name="M8: MICE + SMOTE")
)

ablation_results <- data.frame()
roc_list <- list(); dags_list <- list(); arc_str_list <- list()
cm_list <- list(); fit_list <- list()

for (ex in experiments) {
  tryCatch({
    r <- evaluate_variant(ex$train, ex$val, ex$test, ex$name, boot_R = 200)
    ablation_results <- rbind(ablation_results, r$results)
    roc_list[[ex$name]] <- r$roc; dags_list[[ex$name]] <- r$structure
    arc_str_list[[ex$name]] <- r$arc_strength; cm_list[[ex$name]] <- r$confusion_matrix
    fit_list[[ex$name]] <- r$fit
  }, error=function(e) cat(sprintf("Error in %s: %s\n", ex$name, e$message)))
}

# 95% CI for AUC (bootstrap), then RANK BY AUC (primary), Youden's J (secondary).
master_leaderboard <- ablation_results
master_leaderboard$Test_AUC_95_CI <- vapply(master_leaderboard$Model, function(m){
  ci <- ci.auc(roc_list[[m]], method="bootstrap", boot.n=2000, quiet=TRUE)
  sprintf("[%.4f - %.4f]", ci[1], ci[3]) }, character(1))

master_leaderboard <- master_leaderboard %>%
  relocate(Test_AUC_95_CI, .after = Test_AUC) %>%
  arrange(desc(Test_AUC), desc(Test_Youden_J))     # <-- AUC-primary ranking

print(master_leaderboard)

                                     Model Val_Threshold Test_AUC    Test_AUC_95_CI Test_Sensitivity Test_Specificity
Sensitivity5         M6: MICE + Continuous        0.0442   0.8115 [0.7601 - 0.8566]             0.78           0.6872
Sensitivity2          M3: Continuous Nodes        0.0481   0.8031 [0.7548 - 0.8479]             0.78           0.6944
Sensitivity4 M5: ALL (MICE + Cont + SMOTE)        0.3247   0.7836 [0.7193 - 0.8376]             0.84           0.6481
Sensitivity6        M7: Continuous + SMOTE        0.5985   0.7808 [0.7173 - 0.8422]             0.56           0.7726
Sensitivity1         M2: MICE + Discretize        0.0430   0.7200 [0.6525 - 0.7840]             0.62           0.6893
Sensitivity                   M1: Baseline        0.0550   0.7176 [0.6466 - 0.7860]             0.56           0.7438
Sensitivity7              M8: MICE + SMOTE        0.5290   0.6921 [0.6150 - 0.7636]             0.52           0.7582
Sensitivity3         M4: SMOTE on Baseline        0.1770   0.6841 [0.6034 - 0.7576]             0.66           0.6091
             Test_Youden_J Test_F1
Sensitivity5        0.4672  0.1985
Sensitivity2        0.4744  0.2021
Sensitivity4        0.4881  0.1935
Sensitivity6        0.3326  0.1873
Sensitivity1        0.3093  0.1619
Sensitivity         0.3038  0.1713
Sensitivity7        0.2782  0.1672
Sensitivity3        0.2691  0.1425

write.csv(master_leaderboard, "leaderboard_final.csv", row.names = FALSE)

#DAG structure, CPT and inference table # ============================================================ # SEGMENT 9.5 — MANUSCRIPT ARTIFACTS (DAG, CPTs, examples) # Everything is extracted from the optimal model M6 and written to # CSV/PDF so the manuscript tables/figures use exact, reproducible values. # ============================================================

library(bnlearn); library(dplyr)
BEST <- "M6: MICE + Continuous"
net  <- fit_list[[BEST]]        # fitted hybrid network
dag  <- dags_list[[BEST]]       # structure only
astr <- arc_str_list[[BEST]]    # bootstrap arc strengths
cat("Parents of stroke:", paste(parents(net, "stroke"), collapse=", "), "\n")

Parents of stroke: hypertension, heart_disease, ever_married 

cat("Children of stroke:", paste(children(net, "stroke"), collapse=", "), "\n")

Children of stroke: age, avg_glucose_level 

9.5a DAG — arc list, bootstrap strengths, and figure

# (i) Arc list with bootstrap support (strength) and direction confidence.
arc_tbl <- astr %>%
  filter(strength > 0.50 & direction >= 0.50) %>%
  arrange(desc(strength))
print(arc_tbl)

            from                to strength direction
1   hypertension            stroke    1.000 1.0000000
2  heart_disease            stroke    1.000 1.0000000
3   ever_married         work_type    1.000 0.5300000
4      work_type               bmi    1.000 1.0000000
5      work_type    smoking_status    1.000 0.5600000
6         stroke               age    1.000 1.0000000
7         stroke avg_glucose_level    1.000 1.0000000
8            bmi avg_glucose_level    0.980 0.8010204
9   ever_married avg_glucose_level    0.970 1.0000000
10           age avg_glucose_level    0.960 0.9192708
11           age               bmi    0.945 0.5449735
12  ever_married               age    0.895 1.0000000
13  ever_married     heart_disease    0.875 0.7028571
14  ever_married            stroke    0.875 1.0000000
15     work_type      hypertension    0.670 0.6604478
16        gender     heart_disease    0.530 1.0000000
17     work_type               age    0.520 1.0000000

write.csv(arc_tbl, "M6_arc_strengths.csv", row.names = FALSE)

# (ii) Publication DAG with edge thickness ~ bootstrap strength.
pdf("M6_DAG.pdf", width = 9, height = 7)
strength.plot(dag, astr, main = "M6: MICE + Continuous — consensus DAG",
              shape = "ellipse")
dev.off()

null device 
          1 

print(dag)                                  # full structure summary (nodes, arcs, parents/children)

  Bayesian network learned via Score-based methods

  model:
   [gender][Residence_type][heart_disease|gender][ever_married|heart_disease][work_type|ever_married]
   [hypertension|work_type][smoking_status|work_type][stroke|hypertension:heart_disease:ever_married]
   [age|ever_married:work_type:stroke][bmi|age:work_type][avg_glucose_level|age:ever_married:bmi:stroke]
  nodes:                                 11 
  arcs:                                  17 
    undirected arcs:                     0 
    directed arcs:                       17 
  average markov blanket size:           4.00 
  average neighbourhood size:            3.09 
  average branching factor:              1.55 

  learning algorithm:                    Hill-Climbing 
  score:                                 BIC (cond. Gauss.) 
  penalization coefficient:              4.013901 
  tests used in the learning procedure:  211 
  optimized:                             TRUE 

cat("Model string:\n", modelstring(dag), "\n\n")   # compact one-line DAG

Model string:
 [gender][Residence_type][heart_disease|gender][ever_married|heart_disease][work_type|ever_married][hypertension|work_type][smoking_status|work_type][stroke|hypertension:heart_disease:ever_married][age|ever_married:work_type:stroke][bmi|age:work_type][avg_glucose_level|age:ever_married:bmi:stroke] 

cat("Arc list (from -> to):\n"); print(arcs(dag))  # every directed edge

Arc list (from -> to):
      from            to                 
 [1,] "stroke"        "age"              
 [2,] "heart_disease" "stroke"           
 [3,] "hypertension"  "stroke"           
 [4,] "stroke"        "avg_glucose_level"
 [5,] "work_type"     "age"              
 [6,] "ever_married"  "work_type"        
 [7,] "work_type"     "bmi"              
 [8,] "work_type"     "smoking_status"   
 [9,] "ever_married"  "age"              
[10,] "ever_married"  "avg_glucose_level"
[11,] "work_type"     "hypertension"     
[12,] "age"           "avg_glucose_level"
[13,] "age"           "bmi"              
[14,] "bmi"           "avg_glucose_level"
[15,] "heart_disease" "ever_married"     
[16,] "gender"        "heart_disease"    
[17,] "ever_married"  "stroke"           

cat("\nMarkov blanket of stroke:", paste(mb(dag, "stroke"), collapse = ", "), "\n")

Markov blanket of stroke: age, hypertension, heart_disease, ever_married, work_type, avg_glucose_level, bmi 

cat("Saved M6_DAG.pdf\n")

Saved M6_DAG.pdf

# (iii) Optional interactive DAG (uncomment if visNetwork is installed)
 library(visNetwork)
 nodes <- data.frame(id = nodes(dag), label = nodes(dag))
 edges <- data.frame(from = arcs(dag)[,1], to = arcs(dag)[,2], arrows = "to")
 visNetwork(nodes, edges) %>% visEdges(smooth = FALSE)

9.5b Table 5 — CPT of stroke (discrete parents)

# Stroke is discrete with discrete parents (CGBN forbids continuous parents),
# so its CPT is a clean probability table. We tidy it and flag the key rows.
stroke_cpt <- as.data.frame.table(coef(net$stroke), responseName = "prob")

# Keep P(stroke = 1 | parents) and order by risk.
risk1 <- stroke_cpt %>% filter(stroke == "1") %>%
  mutate(prob_pct = round(100 * prob, 2)) %>%
  arrange(desc(prob)) %>% select(-stroke)
print(risk1)

  hypertension heart_disease ever_married       prob prob_pct
1            1             1           No 0.75000000    75.00
2            0             1           No 0.15384615    15.38
3            1             0           No 0.13333333    13.33
4            1             1          Yes 0.13333333    13.33
5            0             1          Yes 0.12962963    12.96
6            1             0          Yes 0.10878661    10.88
7            0             0          Yes 0.05137615     5.14
8            0             0           No 0.01192843     1.19

write.csv(risk1, "M6_CPT_stroke.csv", row.names = FALSE)

# Auto-extract the worked examples the manuscript quotes:
cat("\n--- Worked CPT examples (exact values for the text) ---\n")

--- Worked CPT examples (exact values for the text) ---

cat(sprintf("Highest-risk parent combination: %.2f%%\n", max(risk1$prob_pct)))

Highest-risk parent combination: 75.00%

cat(sprintf("Lowest-risk parent combination:  %.2f%%\n", min(risk1$prob_pct)))

Lowest-risk parent combination:  1.19%

# Print the full label of the highest-risk cell so you can describe it precisely:
top <- risk1[which.max(risk1$prob_pct), ]
cat("Highest-risk cell:\n"); print(top)

Highest-risk cell:
  hypertension heart_disease ever_married prob prob_pct
1            1             1           No 0.75       75

9.5c Table 6 — continuous children of stroke (Gaussian parameters)

# For each continuous child of stroke, report the local regression coefficients
# and residual SD per discrete-parent configuration (this is the CG parameterisation).
cont_children <- intersect(children(net, "stroke"),
                           c("age", "avg_glucose_level", "bmi"))

for (nd in cont_children) {
  cat("\n==== Continuous node:", nd, "====\n")
  cat("Coefficients (intercept/slopes per parent configuration):\n")
  print(net[[nd]]$coefficients)
  cat("Residual SD per configuration:\n")
  print(net[[nd]]$sd)
  # Tidy export
  co <- as.data.frame(net[[nd]]$coefficients)
  co$term <- rownames(co)
  write.csv(co, sprintf("M6_params_%s_coef.csv", nd), row.names = FALSE)
  write.csv(data.frame(config = names(net[[nd]]$sd), sd = as.numeric(net[[nd]]$sd)),
            sprintf("M6_params_%s_sd.csv", nd), row.names = FALSE)
}

==== Continuous node: age ====
Coefficients (intercept/slopes per parent configuration):
                   0  1        2       3    4  5        6        7        8        9   10 11    12     13 14 15   16
(Intercept) 6.732127 NA 35.90164 52.3141 16.8 NA 27.76749 50.43486 48.77273 59.87949 7.66 NA 69.25 68.375 NA NA 71.2
                  17   18       19
(Intercept) 65.35443 66.4 71.21212
Residual SD per configuration:
        0         1         2         3         4         5         6         7         8         9        10 
 4.450686        NA 15.831093 13.120909  2.573368        NA 13.648362 15.352049 24.172327 15.116798  8.966114 
       11        12        13        14        15        16        17        18        19 
       NA  9.878428  9.500000        NA        NA 10.261037 12.535118 16.652327 10.270582 

==== Continuous node: avg_glucose_level ====
Coefficients (intercept/slopes per parent configuration):
                      0          1          2          3
(Intercept) 89.55712839 37.2230332 13.3129777 -11.671535
age          0.29773939  0.6970657  0.4969107   1.027348
bmi          0.01582724  1.1373560  2.3704633   2.499728
Residual SD per configuration:
       0        1        2        3 
31.89516 46.49507 53.83766 59.86714 

cat("\nSaved per-node coefficient and SD CSVs.\n")

Saved per-node coefficient and SD CSVs.

9.5d Inference & counterfactual examples (Table for the text)

# Exact, reproducible probabilities for the clinical personas + intervention.
set.seed(2026); n_sim <- 1e6   # logic sampling

pA <- cpquery(net, event=(stroke=="1"),
              evidence=(age<=50 & hypertension=="0" & heart_disease=="0"), n=n_sim)
pB <- cpquery(net, event=(stroke=="1"),
              evidence=(age>=65 & hypertension=="1" & heart_disease=="1" & avg_glucose_level>=150), n=n_sim)
pC <- cpquery(net, event=(stroke=="1"),
              evidence=(age>=65 & hypertension=="1" & heart_disease=="1" & avg_glucose_level<=100), n=n_sim)

inference_tbl <- data.frame(
  Scenario = c("A: Low-risk (age<=50, no HTN, no HD)",
               "B: High-risk (age>=65, HTN, HD, glucose>=150)",
               "C: B + glycaemic control (glucose<=100)"),
  Stroke_Probability_pct = round(100 * c(pA, pB, pC), 2))
print(inference_tbl)

                                       Scenario Stroke_Probability_pct
1          A: Low-risk (age<=50, no HTN, no HD)                   0.74
2 B: High-risk (age>=65, HTN, HD, glucose>=150)                  43.18
3       C: B + glycaemic control (glucose<=100)                  31.55

cat(sprintf("\nAbsolute risk reduction (B -> C): %.2f%% -> %.2f%% (delta = %.2f points)\n",
            100*pB, 100*pC, 100*(pB - pC)))

Absolute risk reduction (B -> C): 43.18% -> 31.55% (delta = 11.63 points)

write.csv(inference_tbl, "M6_inference_examples.csv", row.names = FALSE)

# Backward (diagnostic) inference: expected glucose given a stroke occurred.
set.seed(42)
post <- cpdist(net, nodes = "avg_glucose_level",
               evidence = list(stroke = factor("1", levels = levels(train_m6$stroke))),
               method = "lw", n = 1e5)
cat(sprintf("Expected avg glucose among stroke cases (model): %.2f mg/dL\n",
            mean(post$avg_glucose_level, na.rm = TRUE)))

Expected avg glucose among stroke cases (model): 123.59 mg/dL

10 ROC figure

cbPalette <- c("#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7")
# Plot on a DISPLAY COPY so roc_list keeps its full model names (needed in Seg 13).
roc_plot <- roc_list
names(roc_plot) <- gsub(" \\(.*\\)", "", names(roc_plot))
p <- pROC::ggroc(roc_plot, legacy.axes=TRUE, size=1.0) +
  geom_abline(intercept=0, slope=1, colour="darkgray", linetype="dashed") +
  scale_color_manual(values=cbPalette) +
  labs(x="False Positive Rate (1 - Specificity)", y="True Positive Rate (Sensitivity)",
       title="Multi-model ROC comparison") +
  theme_minimal() +
  theme(legend.title=element_blank(), legend.position=c(0.80,0.20))
ggsave("Stroke_Prediction_Multi_ROC.png", p, width=8, height=6, dpi=300)
print(p)

11 Calibration & Brier score

# Brier score = mean((p - y)^2): lower is better-calibrated. Demonstrates the
# claim that SMOTE distorts probabilities (it should INFLATE the Brier score).
brier <- function(m){
  y <- as.numeric(as.character(roc_list[[m]]$response))
  p <- roc_list[[m]]$predictor
  mean((p - y)^2)
}
brier_tbl <- data.frame(Model=names(roc_list),
                        Brier=round(vapply(names(roc_list), brier, numeric(1)),4)) %>%
  arrange(Brier)
print(brier_tbl); write.csv(brier_tbl, "brier_scores.csv", row.names=FALSE)

                                                      Model  Brier
M1: Baseline                                   M1: Baseline 0.0506
M2: MICE + Discretize                 M2: MICE + Discretize 0.0506
M6: MICE + Continuous                 M6: MICE + Continuous 0.0558
M3: Continuous Nodes                   M3: Continuous Nodes 0.0563
M4: SMOTE on Baseline                 M4: SMOTE on Baseline 0.1734
M8: MICE + SMOTE                           M8: MICE + SMOTE 0.1754
M5: ALL (MICE + Cont + SMOTE) M5: ALL (MICE + Cont + SMOTE) 0.1824
M7: Continuous + SMOTE               M7: Continuous + SMOTE 0.1836

# Reliability curve for the four headline arms.
# Tie-robust binning: discrete/SMOTE models produce many identical probabilities,
# so quantile cut-points collapse. We use UNIQUE quantile breaks and fall back to
# binning by distinct probability values when ties are severe.
calib_curve <- function(m, bins=10){
  y <- as.numeric(as.character(roc_list[[m]]$response)); p <- roc_list[[m]]$predictor
  br <- unique(quantile(p, probs=seq(0,1,length.out=bins+1), na.rm=TRUE))
  if (length(br) < 3) {                      # too many ties -> bin by value
    g <- cut(p, breaks=unique(c(-Inf, sort(unique(p)))), include.lowest=TRUE)
  } else {
    g <- cut(p, breaks=br, include.lowest=TRUE)
  }
  out <- data.frame(model=m, pred=tapply(p, g, mean), obs=tapply(y, g, mean))
  na.omit(out)
}
focus <- intersect(c("M1: Baseline","M3: Continuous Nodes","M6: MICE + Continuous",
                     "M5: ALL (MICE + Cont + SMOTE)"), names(roc_list))
cal <- bind_rows(lapply(focus, calib_curve))
ggplot(cal, aes(pred, obs, colour=model)) +
  geom_abline(slope=1, intercept=0, linetype="dashed", colour="grey50") +
  geom_line() + geom_point() +
  labs(x="Mean predicted probability", y="Observed stroke fraction",
       title="Calibration (reliability) curves") + theme_minimal() -> cal_plot
ggsave("calibration_curves.png", cal_plot, width=7, height=6, dpi=300)
print(cal_plot)

12 Opaque-model benchmark: logistic regression & XGBoost

# Same splits, continuous features. Quantifies how close the interpretable BN
# (AUC ~0.81) comes to standard black-box models — the paper's motivating tension.
suppressPackageStartupMessages({library(xgboost)})

bench_tr <- train_m3; bench_va <- val_m3; bench_te <- test_m3   # continuous, no SMOTE

# --- Logistic regression
lr <- glm(stroke ~ ., data=bench_tr, family=binomial)
lr_te <- predict(lr, bench_te, type="response")
lr_auc <- as.numeric(auc(roc(bench_te$stroke, lr_te, levels=c("0","1"), quiet=TRUE)))

# --- XGBoost (stable xgb.train API; compatible with xgboost 1.x-3.x)
# nrounds is tuned by EARLY STOPPING on the validation fold so the comparison is
# fair (a fixed large nrounds overfits on ~249 events and understates XGBoost).
mm_tr <- model.matrix(stroke ~ . -1, data=bench_tr)
mm_va <- model.matrix(stroke ~ . -1, data=bench_va)
mm_te <- model.matrix(stroke ~ . -1, data=bench_te)
feat  <- Reduce(intersect, list(colnames(mm_tr), colnames(mm_va), colnames(mm_te)))
dtr <- xgb.DMatrix(mm_tr[, feat, drop=FALSE], label=as.numeric(as.character(bench_tr$stroke)))
dval<- xgb.DMatrix(mm_va[, feat, drop=FALSE], label=as.numeric(as.character(bench_va$stroke)))
dte <- xgb.DMatrix(mm_te[, feat, drop=FALSE], label=as.numeric(as.character(bench_te$stroke)))
params <- list(objective="binary:logistic", eval_metric="auc",
               max_depth=3, learning_rate=0.05, subsample=0.8, colsample_bytree=0.8,
               scale_pos_weight=sum(bench_tr$stroke=="0")/sum(bench_tr$stroke=="1"))
set.seed(GLOBAL_SEED)
xgb <- xgb.train(params=params, data=dtr, nrounds=1000,
                 watchlist=list(val=dval), early_stopping_rounds=25, verbose=0)
xgb_te <- predict(xgb, dte)
xgb_auc <- as.numeric(auc(roc(bench_te$stroke, xgb_te, levels=c("0","1"), quiet=TRUE)))

bench_tbl <- data.frame(
  Model=c("Bayesian network (M6)","Logistic regression","XGBoost"),
  Test_AUC=round(c(master_leaderboard$Test_AUC[master_leaderboard$Model=="M6: MICE + Continuous"],
                   lr_auc, xgb_auc),4))
print(bench_tbl); write.csv(bench_tbl, "benchmark_opaque.csv", row.names=FALSE)

                  Model Test_AUC
1 Bayesian network (M6)   0.8115
2   Logistic regression   0.8188
3               XGBoost   0.8222

13 Hypothesis testing — DeLong (AUC) + McNemar (sensitivity)

# Thresholds are read from the leaderboard (no hard-coding).
thr <- function(m) master_leaderboard$Val_Threshold[master_leaderboard$Model==m]

# McNemar on the actual-stroke cohort: did model A catch cases B missed?
mcnemar_sens <- function(A, B){
  y  <- roc_list[[A]]$response
  pa <- ifelse(roc_list[[A]]$predictor >= thr(A), 1, 0)
  pb <- ifelse(roc_list[[B]]$predictor >= thr(B), 1, 0)
  idx <- which(y == 1)
  tab <- table(A=factor(pa[idx],0:1), B=factor(pb[idx],0:1))
  cat(sprintf("\nMcNemar (sensitivity) %s vs %s:\n", A, B)); print(mcnemar.test(tab))
}

cat("=== H1 Topology: continuous vs discrete ===\n")

=== H1 Topology: continuous vs discrete ===

print(roc.test(roc_list[["M3: Continuous Nodes"]], roc_list[["M1: Baseline"]], method="delong"))

    DeLong's test for two correlated ROC curves

data:  roc_list[["M3: Continuous Nodes"]] and roc_list[["M1: Baseline"]]
Z = 2.2989, p-value = 0.02151
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 0.01261342 0.15849769
sample estimates:
AUC of roc1 AUC of roc2 
  0.8031070   0.7175514 

print(roc.test(roc_list[["M6: MICE + Continuous"]], roc_list[["M2: MICE + Discretize"]], method="delong"))

    DeLong's test for two correlated ROC curves

data:  roc_list[["M6: MICE + Continuous"]] and roc_list[["M2: MICE + Discretize"]]
Z = 2.546, p-value = 0.0109
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 0.02106229 0.16194183
sample estimates:
AUC of roc1 AUC of roc2 
  0.8114815   0.7199794 

cat("\n=== H2 Imputation: MICE vs median ===\n")

=== H2 Imputation: MICE vs median ===

print(roc.test(roc_list[["M2: MICE + Discretize"]], roc_list[["M1: Baseline"]], method="delong"))

    DeLong's test for two correlated ROC curves

data:  roc_list[["M2: MICE + Discretize"]] and roc_list[["M1: Baseline"]]
Z = 0.37426, p-value = 0.7082
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 -0.01028730  0.01514327
sample estimates:
AUC of roc1 AUC of roc2 
  0.7199794   0.7175514 

print(roc.test(roc_list[["M6: MICE + Continuous"]], roc_list[["M3: Continuous Nodes"]], method="delong"))

    DeLong's test for two correlated ROC curves

data:  roc_list[["M6: MICE + Continuous"]] and roc_list[["M3: Continuous Nodes"]]
Z = 1.3128, p-value = 0.1893
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 -0.004128663  0.020877634
sample estimates:
AUC of roc1 AUC of roc2 
  0.8114815   0.8031070 

cat("\n=== H3 Balancing: SMOTE vs unadjusted (AUC + sensitivity) ===\n")

=== H3 Balancing: SMOTE vs unadjusted (AUC + sensitivity) ===

# Key point: SMOTE does NOT improve AUC (DeLong), even where it raises sensitivity.
M5 <- "M5: ALL (MICE + Cont + SMOTE)"   # full name, matches roc_list / leaderboard
print(roc.test(roc_list[["M4: SMOTE on Baseline"]], roc_list[["M1: Baseline"]], method="delong"))

    DeLong's test for two correlated ROC curves

data:  roc_list[["M4: SMOTE on Baseline"]] and roc_list[["M1: Baseline"]]
Z = -1.8832, p-value = 0.05968
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 -0.068278205  0.001364625
sample estimates:
AUC of roc1 AUC of roc2 
  0.6840947   0.7175514 

print(roc.test(roc_list[[M5]], roc_list[["M6: MICE + Continuous"]], method="delong"))

    DeLong's test for two correlated ROC curves

data:  roc_list[[M5]] and roc_list[["M6: MICE + Continuous"]]
Z = -1.07, p-value = 0.2846
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 -0.07889048  0.02317031
sample estimates:
AUC of roc1 AUC of roc2 
  0.7836214   0.8114815 

mcnemar_sens(M5, "M6: MICE + Continuous")

McNemar (sensitivity) M5: ALL (MICE + Cont + SMOTE) vs M6: MICE + Continuous:

    McNemar's Chi-squared test with continuity correction

data:  tab
McNemar's chi-squared = 0.57143, df = 1, p-value = 0.4497

cat("\n=== H4 Synergy: best continuous vs absolute baseline ===\n")

=== H4 Synergy: best continuous vs absolute baseline ===

print(roc.test(roc_list[["M6: MICE + Continuous"]], roc_list[["M1: Baseline"]], method="delong"))

    DeLong's test for two correlated ROC curves

data:  roc_list[["M6: MICE + Continuous"]] and roc_list[["M1: Baseline"]]
Z = 2.5075, p-value = 0.01216
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 0.02051173 0.16734835
sample estimates:
AUC of roc1 AUC of roc2 
  0.8114815   0.7175514 

#13b Hypothesis test and p-value re-run

library(pROC)

# Validation-tuned threshold for each model (read from the leaderboard).
thr <- function(m) master_leaderboard$Val_Threshold[master_leaderboard$Model == m]

# McNemar on the actual-stroke cohort: does model A catch cases B misses?
mcnemar_sens <- function(A, B){
  y  <- roc_list[[A]]$response
  pa <- ifelse(roc_list[[A]]$predictor >= thr(A), 1, 0)
  pb <- ifelse(roc_list[[B]]$predictor >= thr(B), 1, 0)
  idx <- which(y == 1)
  tab <- table(A = factor(pa[idx], 0:1), B = factor(pb[idx], 0:1))
  p <- suppressWarnings(mcnemar.test(tab)$p.value)
  cat(sprintf("McNemar sens  %-30s vs %-24s : p = %.4f\n", A, B, p))
  data.frame(test = "McNemar(sens)", A = A, B = B, stat = NA, p = round(p, 4))
}

# DeLong AUC comparison.
delong <- function(A, B){
  t <- roc.test(roc_list[[A]], roc_list[[B]], method = "delong")
  cat(sprintf("DeLong AUC    %-30s vs %-24s : Z = %.3f, p = %.4f\n",
              A, B, as.numeric(t$statistic), t$p.value))
  data.frame(test = "DeLong(AUC)", A = A, B = B,
             stat = round(as.numeric(t$statistic), 3), p = round(t$p.value, 4))
}

M5 <- "M5: ALL (MICE + Cont + SMOTE)"
res <- list()

cat("\n=== H1  topology, minority sensitivity (McNemar) ===\n")

=== H1  topology, minority sensitivity (McNemar) ===

res[[1]] <- mcnemar_sens("M3: Continuous Nodes",  "M1: Baseline")

McNemar sens  M3: Continuous Nodes           vs M1: Baseline             : p = 0.0098

res[[2]] <- mcnemar_sens("M6: MICE + Continuous", "M2: MICE + Discretize")

McNemar sens  M6: MICE + Continuous          vs M2: MICE + Discretize    : p = 0.0990

cat("\n=== H3  SMOTE, minority sensitivity (McNemar) ===\n")

=== H3  SMOTE, minority sensitivity (McNemar) ===

res[[3]] <- mcnemar_sens("M4: SMOTE on Baseline", "M1: Baseline")

McNemar sens  M4: SMOTE on Baseline          vs M1: Baseline             : p = 0.1824

res[[4]] <- mcnemar_sens("M8: MICE + SMOTE",      "M2: MICE + Discretize")

McNemar sens  M8: MICE + SMOTE               vs M2: MICE + Discretize    : p = 0.2278

res[[5]] <- mcnemar_sens(M5,                       "M6: MICE + Continuous")

McNemar sens  M5: ALL (MICE + Cont + SMOTE)  vs M6: MICE + Continuous    : p = 0.4497

cat("\n=== H4  synergy vs absolute baseline (DeLong AUC) ===\n")

=== H4  synergy vs absolute baseline (DeLong AUC) ===

res[[6]] <- delong(M5,                       "M1: Baseline")

DeLong AUC    M5: ALL (MICE + Cont + SMOTE)  vs M1: Baseline             : Z = 1.767, p = 0.0772

res[[7]] <- delong("M6: MICE + Continuous",  "M1: Baseline")

DeLong AUC    M6: MICE + Continuous          vs M1: Baseline             : Z = 2.508, p = 0.0122

pval_table <- do.call(rbind, res)
print(pval_table)

           test                             A                     B  stat      p
1 McNemar(sens)          M3: Continuous Nodes          M1: Baseline    NA 0.0098
2 McNemar(sens)         M6: MICE + Continuous M2: MICE + Discretize    NA 0.0990
3 McNemar(sens)         M4: SMOTE on Baseline          M1: Baseline    NA 0.1824
4 McNemar(sens)              M8: MICE + SMOTE M2: MICE + Discretize    NA 0.2278
5 McNemar(sens) M5: ALL (MICE + Cont + SMOTE) M6: MICE + Continuous    NA 0.4497
6   DeLong(AUC) M5: ALL (MICE + Cont + SMOTE)          M1: Baseline 1.767 0.0772
7   DeLong(AUC)         M6: MICE + Continuous          M1: Baseline 2.508 0.0122

write.csv(pval_table, "hypothesis_pvalues.csv", row.names = FALSE)
cat("\nSaved hypothesis_pvalues.csv\n")

Saved hypothesis_pvalues.csv

14 Robustness — 30 repeated stratified splits

# A single split is fragile on 249 events. Repeat the WHOLE pipeline over many
# seeds; report mean +/- SD AUC, and a PAIRED test of continuous vs discrete on
# the per-seed AUC differences. This is far stronger than one DeLong p-value.
# NOTE: continuous arms use per-patient likelihood weighting; n is lowered to 500
# here for runtime. Raise N_REPEATS / n for the final camera-ready run.
N_REPEATS <- 30

eval_light <- function(tr, va, te, disc){
  wl <- if (disc) wl_discrete else wl_hybrid; bl <- if (disc) bl_discrete else bl_hybrid
  dag <- hc(tr, whitelist=wl, blacklist=bl)
  fit <- if (disc) bn.fit(dag, tr, method="bayes", iss=10) else bn.fit(dag, tr)
  gp <- function(d){
    if (disc){ pp <- predict(fit, node="stroke", data=d, method="bayes-lw", prob=TRUE)
               return(attr(pp,"prob")["1",]) }
    ec <- setdiff(names(d),"stroke")
    vapply(seq_len(nrow(d)), function(i){
      p <- suppressWarnings(cpquery(fit, event=(stroke=="1"),
             evidence=as.list(d[i,ec]), method="lw", n=500)); if (is.na(p)) 0 else p}, numeric(1))
  }
  as.numeric(auc(roc(te$stroke, gp(te), levels=c("0","1"), quiet=TRUE)))
}

rob <- data.frame()
for (s in 1:N_REPEATS) {
  sp <- split_stroke(stroke_data, seed = 1000 + s)
  sp$train$split<-NULL; # safety
  m  <- median(as.numeric(as.character(sp$train$bmi)), na.rm=TRUE)
  # MICE (leakage-safe) for this split
  cb <- bind_rows(mutate(sp$train,split="train"), mutate(sp$val,split="val"),
                  mutate(sp$test,split="test")) %>% select(-any_of("id"))
  mm <- suppressWarnings(mice(cb, m=1, ignore=cb$split!="train", maxit=5,
                              seed=1000+s, printFlag=FALSE))
  dn <- complete(mm,1)
  tr2<-dn%>%filter(split=="train")%>%select(-split)
  te2<-dn%>%filter(split=="test") %>%select(-split)
  safe <- function(expr) tryCatch(expr, error=function(e){ cat("  (skipped:",e$message,")\n"); NA_real_ })
  rob <- rbind(rob, data.frame(seed=s,
    M1=safe(eval_light(impute_baseline(sp$train,m), NULL, impute_baseline(sp$test,m), TRUE)),
    M2=safe(eval_light(impute_baseline(tr2,m),      NULL, impute_baseline(te2,m),      TRUE)),
    M3=safe(eval_light(prep_continuous(sp$train,m), NULL, prep_continuous(sp$test,m),  FALSE)),
    M6=safe(eval_light(prep_continuous(tr2,m),      NULL, prep_continuous(te2,m),       FALSE))))
  cat(sprintf("robust seed %d/%d done\n", s, N_REPEATS))
}

robust seed 1/30 done
robust seed 2/30 done
robust seed 3/30 done
robust seed 4/30 done
robust seed 5/30 done
robust seed 6/30 done
robust seed 7/30 done
robust seed 8/30 done
robust seed 9/30 done
robust seed 10/30 done
robust seed 11/30 done
robust seed 12/30 done
  (skipped: 'Residence_type' has different number of levels in the node and in the data. )
  (skipped: 'Residence_type' has different number of levels in the node and in the data. )
robust seed 13/30 done
robust seed 14/30 done
robust seed 15/30 done
robust seed 16/30 done
robust seed 17/30 done
robust seed 18/30 done
robust seed 19/30 done
robust seed 20/30 done
robust seed 21/30 done
robust seed 22/30 done
robust seed 23/30 done
robust seed 24/30 done
robust seed 25/30 done
robust seed 26/30 done
robust seed 27/30 done
robust seed 28/30 done
robust seed 29/30 done
robust seed 30/30 done

rob_summary <- data.frame(
  Arm = c("M1: Baseline","M2: MICE+Disc","M3: Continuous","M6: MICE+Continuous"),
  AUC_mean = round(c(mean(rob$M1,na.rm=TRUE), mean(rob$M2,na.rm=TRUE), mean(rob$M3,na.rm=TRUE), mean(rob$M6,na.rm=TRUE)),4),
  AUC_sd   = round(c(sd(rob$M1,na.rm=TRUE),   sd(rob$M2,na.rm=TRUE),   sd(rob$M3,na.rm=TRUE),   sd(rob$M6,na.rm=TRUE)),4))
print(rob_summary); write.csv(rob_summary, "robustness_summary.csv", row.names=FALSE)

                  Arm AUC_mean AUC_sd
1        M1: Baseline   0.7696 0.0392
2       M2: MICE+Disc   0.7693 0.0402
3      M3: Continuous   0.8209 0.0284
4 M6: MICE+Continuous   0.8202 0.0283

cat("\nPaired tests (continuous vs discrete) across", N_REPEATS, "seeds:\n")

Paired tests (continuous vs discrete) across 30 seeds:

cat("H1 primary  M3 vs M1:\n"); print(t.test(rob$M3, rob$M1, paired=TRUE)); print(wilcox.test(rob$M3, rob$M1, paired=TRUE))

H1 primary  M3 vs M1:

    Paired t-test

data:  rob$M3 and rob$M1
t = 8.1395, df = 28, p-value = 7.333e-09
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.03899015 0.06521447
sample estimates:
mean difference 
     0.05210231 


    Wilcoxon signed rank exact test

data:  rob$M3 and rob$M1
V = 429, p-value = 5.215e-08
alternative hypothesis: true location shift is not equal to 0

cat("H1 support  M6 vs M2:\n"); print(t.test(rob$M6, rob$M2, paired=TRUE)); print(wilcox.test(rob$M6, rob$M2, paired=TRUE))

H1 support  M6 vs M2:

    Paired t-test

data:  rob$M6 and rob$M2
t = 8.1929, df = 28, p-value = 6.434e-09
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.03875949 0.06460222
sample estimates:
mean difference 
     0.05168086 


    Wilcoxon signed rank exact test

data:  rob$M6 and rob$M2
V = 426, p-value = 1.229e-07
alternative hypothesis: true location shift is not equal to 0

write.csv(rob, "robustness_raw_auc.csv", row.names=FALSE)

15 Clinical inference & what-if scenarios

fitted_net <- fit_list[["M6: MICE + Continuous"]]
set.seed(2026); n_sim <- 1e6   # logic sampling; seed fixes the reported values

# Evidence must be written INLINE: cpquery uses non-standard evaluation and will
# not accept a quoted expression passed through a variable.
pA <- cpquery(fitted_net, event=(stroke=="1"),
              evidence=(age<=50 & hypertension=="0" & heart_disease=="0"), n=n_sim)
pB <- cpquery(fitted_net, event=(stroke=="1"),
              evidence=(age>=65 & hypertension=="1" & heart_disease=="1" & avg_glucose_level>=150), n=n_sim)
pC <- cpquery(fitted_net, event=(stroke=="1"),
              evidence=(age>=65 & hypertension=="1" & heart_disease=="1" & avg_glucose_level<=100), n=n_sim)
cat(sprintf("Persona A (low risk):           %.2f%%\n", 100*pA))

Persona A (low risk):           0.74%

cat(sprintf("Persona B (high risk):          %.2f%%\n", 100*pB))

Persona B (high risk):          43.18%

cat(sprintf("Persona C (B + glucose managed): %.2f%%\n", 100*pC))

Persona C (B + glucose managed): 31.55%

cat(sprintf("Absolute risk reduction:        %.2f%% -> %.2f%%\n", 100*pB, 100*pC))

Absolute risk reduction:        43.18% -> 31.55%

16 Session info (reproducibility record)

writeLines(capture.output(sessionInfo()), "sessionInfo.txt")
sessionInfo()

R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_Canada.utf8  LC_CTYPE=English_Canada.utf8    LC_MONETARY=English_Canada.utf8
[4] LC_NUMERIC=C                    LC_TIME=English_Canada.utf8    

time zone: America/Toronto
tzcode source: internal

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] xgboost_3.2.1.1 vcd_1.4-13      themis_1.0.3    recipes_1.3.1   mice_3.19.0     pROC_1.19.0.1   caret_7.0-1    
 [8] lattice_0.22-7  bnlearn_5.1     lubridate_1.9.4 forcats_1.0.1   stringr_1.6.0   dplyr_1.1.4     purrr_1.2.0    
[15] readr_2.1.6     tidyr_1.3.1     tibble_3.3.0    ggplot2_4.0.1   tidyverse_2.0.0

loaded via a namespace (and not attached):
  [1] Rdpack_2.6.4         rlang_1.1.6          magrittr_2.0.4       otel_0.2.0           e1071_1.7-17        
  [6] compiler_4.5.2       systemfonts_1.3.1    vctrs_0.6.5          reshape2_1.4.5       pkgconfig_2.0.3     
 [11] shape_1.4.6.1        fastmap_1.2.0        backports_1.5.0      labeling_0.4.3       utf8_1.2.6          
 [16] rmarkdown_2.30       prodlim_2025.04.28   effsize_0.8.1        tzdb_0.5.0           pracma_2.4.6        
 [21] graph_1.88.1         nloptr_2.2.1         ragg_1.5.0           xfun_0.57            glmnet_4.1-10       
 [26] jomo_2.7-6           cachem_1.1.0         jsonlite_2.0.0       pan_1.9              broom_1.0.11        
 [31] parallel_4.5.2       R6_2.6.1             bslib_0.9.0          stringi_1.8.7        RColorBrewer_1.1-3  
 [36] parallelly_1.46.1    car_3.1-5            boot_1.3-32          rpart_4.1.24         lmtest_0.9-40       
 [41] jquerylib_0.1.4      Rcpp_1.1.0           iterators_1.0.14     knitr_1.51           future.apply_1.20.1 
 [46] zoo_1.8-15           Matrix_1.7-4         splines_4.5.2        nnet_7.3-20          timechange_0.3.0    
 [51] tidyselect_1.2.1     rstudioapi_0.17.1    abind_1.4-8          yaml_2.3.12          timeDate_4051.111   
 [56] codetools_0.2-20     listenv_0.10.0       plyr_1.8.9           withr_3.0.2          S7_0.2.1            
 [61] evaluate_1.0.5       future_1.69.0        survival_3.8-3       proxy_0.4-29         pillar_1.11.1       
 [66] ggpubr_0.6.3         carData_3.0-6        rsconnect_1.7.0      foreach_1.5.2        stats4_4.5.2        
 [71] reformulas_0.4.3.1   generics_0.1.4       hms_1.1.4            scales_1.4.0         minqa_1.2.8         
 [76] globals_0.18.0       class_7.3-23         glue_1.8.0           ROSE_0.0-4           tools_4.5.2         
 [81] data.table_1.17.8    lme4_1.1-38          ModelMetrics_1.2.2.2 gower_1.0.2          ggsignif_0.6.4      
 [86] rbibutils_2.4        colorspace_2.1-2     ipred_0.9-15         nlme_3.1-168         Formula_1.2-5       
 [91] cli_3.6.5            textshaping_1.0.4    lava_1.8.2           Rgraphviz_2.54.0     gtable_0.3.6        
 [96] rstatix_0.7.3        sass_0.4.10          digest_0.6.39        BiocGenerics_0.56.0  farver_2.1.2        
[101] htmltools_0.5.9      lifecycle_1.0.5      hardhat_1.4.2        mitml_0.4-5          MASS_7.3-65