Hybrid Bayesian Network for Interpretable Stroke Prediction

cite: “P. Gc, R.-J. Singh, and R. Grewal, “Explainable AI for Stroke Prediction: A Hybrid Bayesian Network Approach,” Jul.09, 2026, Preprints: 2026070625. doi: 10.20944/preprints202607.0625.v1.”

updated on: July 28, 2026

00 Setup & reproducibility

library(tidyverse); library(bnlearn); library(caret); library(pROC)
library(mice); library(recipes); library(themis)
library(lme4); library(lmerTest)   # factorial mixed model
library(xgboost)                   # opaque benchmark
library(vcd)                       # training-only association screen

GLOBAL_SEED <- 42
set.seed(GLOBAL_SEED)

# Inference is LIKELIHOOD WEIGHTING for every arm.

INFERENCE_METHOD <- "Monte Carlo likelihood weighting"
N_PARTICLES_PATIENT <- 1000      # per-patient posterior  (matches Methods)
N_PARTICLES_QUERY   <- 1e6       # population-level conditional-risk queries

01 Load & clean data

stroke_data <- read.csv("healthcare-dataset-stroke-data.csv")
stroke_data <- stroke_data %>%
  filter(gender != "Other") %>% droplevels() %>%
  mutate(bmi = as.numeric(as.character(bmi)),
         across(c(stroke,hypertension,heart_disease,gender,ever_married,
                  work_type,Residence_type,smoking_status), as.factor))

PI_POP <- mean(stroke_data$stroke == "1")
N_TOTAL <- nrow(stroke_data)
cat(sprintf("N = %d | stroke PREVALENCE = %.2f%% | missing BMI = %d (%.1f%%)\n",
            N_TOTAL, 100*PI_POP, sum(is.na(stroke_data$bmi)),
            100*mean(is.na(stroke_data$bmi))))

N = 5109 | stroke PREVALENCE = 4.87% | missing BMI = 201 (3.9%)

02 Exploratory data analysis (figures)

# Set eval=TRUE to regenerate the feature-distribution figures.

03 Stratified split, absolute sizes, level harmonisation

split_stroke <- function(data, seed = GLOBAL_SEED, p_train = .6, p_val = .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)),]))
}

# Guarantees every resample is analysable. 
harmonise_levels <- function(tr, va, te) {
  for (v in names(Filter(is.factor, tr))) {
    lv <- Reduce(union, list(levels(factor(tr[[v]])), levels(factor(va[[v]])),
                             levels(factor(te[[v]]))))
    tr[[v]] <- factor(as.character(tr[[v]]), levels = lv)
    va[[v]] <- factor(as.character(va[[v]]), levels = lv)
    te[[v]] <- factor(as.character(te[[v]]), levels = lv)
  }
  list(train = tr, val = va, test = te)
}

sp <- split_stroke(stroke_data); sp <- harmonise_levels(sp$train, sp$val, sp$test)
trainset <- sp$train; valset <- sp$val; testset <- sp$test

# Partition data
split_tbl <- tibble(
  Partition = c("Training","Validation","Test"),
  n      = c(nrow(trainset), nrow(valset), nrow(testset)),
  Events = c(sum(trainset$stroke=="1"), sum(valset$stroke=="1"), sum(testset$stroke=="1"))) %>%
  mutate(Prevalence = sprintf("%.2f%%", 100*Events/n))
print(split_tbl); write.csv(split_tbl, "v3_split_sizes.csv", row.names=FALSE)

# A tibble: 3 × 4
  Partition      n Events Prevalence
  <chr>      <int>  <int> <chr>     
1 Training    3065    149 4.86%     
2 Validation  1022     50 4.89%     
3 Test        1022     50 4.89%     

04 Preprocessing helpers (incl. Borderline-SMOTE)

# --- (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="stroke") {
  df <- droplevels(as.data.frame(df))
  rec <- recipe(formula(paste(target,"~ .")), data=df) %>%
    step_smotenc(all_outcomes(), over_ratio=1, seed=GLOBAL_SEED) %>% prep()
  droplevels(as.data.frame(juice(rec) %>%
    mutate(across(where(is.numeric), as.numeric),
           across(where(is.character), as.factor))))
}

05 Build the 8 ablation study arms (M1-M8)

# outcome retained for TRAINING imputation, excluded for VAL/TEST
# single imputation by chained equations (SICE)
M_IMP <- 1
impute_sice <- function(trainset, valset, testset, seed = GLOBAL_SEED) {
  trainset$.split <- "train"; valset$.split <- "val"; testset$.split <- "test"
  comb <- bind_rows(trainset, valset, testset) %>% select(-any_of("id"))
  dat  <- comb %>% select(-.split); ig <- comb$.split != "train"
  pm_tr <- make.predictorMatrix(dat); pm_te <- pm_tr; pm_te[,"stroke"] <- 0
  a <- complete(suppressWarnings(mice(dat, m=M_IMP, maxit=5, predictorMatrix=pm_tr,
                 ignore=ig, seed=seed,   printFlag=FALSE)),1)
  b <- complete(suppressWarnings(mice(dat, m=M_IMP, maxit=5, predictorMatrix=pm_te,
                 ignore=ig, seed=seed+1, printFlag=FALSE)),1)
  list(train=a[comb$.split=="train",,drop=FALSE],
       val  =b[comb$.split=="val",  ,drop=FALSE],
       test =b[comb$.split=="test", ,drop=FALSE])
}
sice <- impute_sice(trainset, valset, testset)
med  <- median(as.numeric(as.character(trainset$bmi)), na.rm=TRUE)

train_m1 <- impute_baseline(trainset, med); val_m1 <- impute_baseline(valset, med); test_m1 <- impute_baseline(testset, med)
train_m2 <- impute_baseline(sice$train, med); val_m2 <- impute_baseline(sice$val, med); test_m2 <- impute_baseline(sice$test, med)
train_m3 <- prep_continuous(trainset, med); val_m3 <- prep_continuous(valset, med); test_m3 <- prep_continuous(testset, med)
train_m4 <- apply_smotenc(train_m1); val_m4 <- val_m1; test_m4 <- test_m1
train_m5_base <- prep_continuous(sice$train, med)
train_m5 <- apply_smotenc(train_m5_base)
val_m5 <- prep_continuous(sice$val, med); test_m5 <- prep_continuous(sice$test, med)
train_m6 <- train_m5_base; val_m6 <- val_m5; test_m6 <- test_m5
train_m7 <- apply_smotenc(train_m3); val_m7 <- val_m3; test_m7 <- test_m3
# [R1.2c] M8 now uses SICE val/test (v1 wrongly used median-imputed sets)
train_m8 <- apply_smotenc(train_m2); val_m8 <- val_m2; test_m8 <- test_m2
cat("Eight arms (2x2x2) built with leakage-safe imputation.\n")

Eight arms (2x2x2) built with leakage-safe imputation.

06 Association screen on training data, supplementary use only

assoc_tbl <- bind_rows(
  map_dfr(c("gender","hypertension","heart_disease","ever_married","work_type",
            "Residence_type","smoking_status"),
    ~data.frame(Variable=.x, Type="categorical",
      Assoc=round(assocstats(table(train_m3[[.x]], train_m3$stroke))$cramer,4))),
  map_dfr(c("age","avg_glucose_level","bmi"),
    ~data.frame(Variable=.x, Type="continuous",
      Assoc=round(cor(as.numeric(train_m3[[.x]]),
        as.numeric(as.character(train_m3$stroke)), use="complete.obs"),4)))
) %>% arrange(desc(abs(Assoc)))
print(assoc_tbl); write.csv(assoc_tbl,"v3_supp_assoc_TRAIN_ONLY.csv",row.names=FALSE)

            Variable        Type  Assoc
1                age  continuous 0.2410
2  avg_glucose_level  continuous 0.1318
3       hypertension categorical 0.1132
4      heart_disease categorical 0.1071
5       ever_married categorical 0.0964
6          work_type categorical 0.0926
7     smoking_status categorical 0.0750
8                bmi  continuous 0.0356
9     Residence_type categorical 0.0250
10            gender categorical 0.0008

07 Structural constraints — specified a priori from the literature

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"))

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")
) %>% as.matrix(); colnames(bl_discrete) <- c("from","to")

# continuous nodes cannot parent discrete nodes -> 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")
) %>% as.matrix(); colnames(bl_hybrid) <- c("from","to")
cat("Whitelist and blacklist create based on literature. \n")

Whitelist and blacklist create based on literature. 

08 Core evaluation function

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)

  # [R1.5b] likelihood weighting for BOTH arm types
  get_probs <- function(td, fit=fitted, disc=is_disc) {
    if (disc) { pp <- predict(fit, node="stroke", data=td, method="bayes-lw", prob=TRUE)
                return(attr(pp,"prob")["1",]) }
    ec <- setdiff(names(td),"stroke")
    vapply(seq_len(nrow(td)), function(i){
      p <- suppressWarnings(cpquery(fit, event=(stroke=="1"),
             evidence=as.list(td[i,ec]), method="lw", n=N_PARTICLES_PATIENT))
      if (is.na(p)) 0 else p}, numeric(1))
  }

  vp <- get_probs(val_df); roc_v <- roc(val_df$stroke, vp, levels=c("0","1"), quiet=TRUE)
  th <- coords(roc_v,"best",ret="threshold",best.method="youden")$threshold[1]
  if (is.na(th)) th <- 0.5
  tp <- get_probs(test_df); roc_t <- roc(test_df$stroke, tp, levels=c("0","1"), quiet=TRUE)
  cm <- confusionMatrix(factor(ifelse(tp>=th,"1","0"),levels=c("0","1")),
                        test_df$stroke, positive="1")

  list(results=data.frame(Model=model_name, Val_Threshold=round(th,4),
         Val_AUC=round(as.numeric(auc(roc_v)),4),
         Test_AUC=round(as.numeric(auc(roc_t)),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_t, roc_val=roc_v, fit=fitted, structure=dag, arc_strength=arc_str,
       confusion_matrix=cm, is_disc=is_disc, wl=wl, bl=bl, train=train_df, test=test_df)
}

09 Main result — train arms, then select on VALIDATION

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

ablation <- data.frame(); roc_list <- list(); dags <- list()
arcs_l <- list(); cms <- list(); fits <- list(); runs <- list()
for (ex in experiments) {
  tryCatch({
    r <- evaluate_variant(ex[[1]],ex[[2]],ex[[3]],ex[[4]],boot_R=200)
    ablation <- rbind(ablation, r$results); roc_list[[ex[[4]]]] <- r$roc
    dags[[ex[[4]]]] <- r$structure; arcs_l[[ex[[4]]]] <- r$arc_strength
    cms[[ex[[4]]]] <- r$confusion_matrix; fits[[ex[[4]]]] <- r$fit; runs[[ex[[4]]]] <- r
  }, error=function(e) cat(sprintf("Error %s: %s\n", ex[[4]], e$message)))
}

leaderboard <- ablation
leaderboard$Test_AUC_95_CI <- vapply(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))

# selection on VALIDATION only, with the 0.005 parsimony rule
vr <- leaderboard %>% arrange(desc(Val_AUC))
near <- vr %>% filter(Val_AUC >= max(Val_AUC) - 0.005)
simple_rank <- c("M1","M3","M2","M6","M4","M7","M8","M5")
near$simp <- match(sub(":.*","",near$Model), simple_rank)
PRIMARY <- near$Model[which.min(near$simp)]
cat("\nBest validation AUC:", vr$Model[1], sprintf("(%.4f)", vr$Val_AUC[1]),
    "\nPRIMARY after 0.005 parsimony rule:", PRIMARY,
    "\nIts unbiased TEST AUC:", leaderboard$Test_AUC[leaderboard$Model==PRIMARY], "\n")

Best validation AUC: M6: SICE + Continuous (0.8683) 
PRIMARY after 0.005 parsimony rule: M3: Continuous Nodes 
Its unbiased TEST AUC: 0.8109 

leaderboard <- leaderboard %>% relocate(Test_AUC_95_CI,.after=Test_AUC) %>%
  arrange(desc(Test_AUC), desc(Test_Youden_J))
print(leaderboard); write.csv(leaderboard,"v3_leaderboard.csv",row.names=FALSE)

                                     Model Val_Threshold Val_AUC Test_AUC    Test_AUC_95_CI Test_Sensitivity
Sensitivity2          M3: Continuous Nodes        0.0857  0.8642   0.8109 [0.7612 - 0.8533]             0.62
Sensitivity5         M6: SICE + Continuous        0.0951  0.8683   0.8094 [0.7615 - 0.8522]             0.58
Sensitivity6        M7: Continuous + SMOTE        0.4303  0.7602   0.7456 [0.6672 - 0.8126]             0.72
Sensitivity4 M5: ALL (SICE + Cont + SMOTE)        0.3022  0.7500   0.7381 [0.6688 - 0.8048]             0.76
Sensitivity                   M1: Baseline        0.0350  0.7710   0.7235 [0.6542 - 0.7907]             0.68
Sensitivity1         M2: SICE + Discretize        0.0330  0.7774   0.7108 [0.6376 - 0.7805]             0.68
Sensitivity3         M4: SMOTE on Baseline        0.5530  0.7509   0.6903 [0.6128 - 0.7630]             0.50
Sensitivity7              M8: SICE + SMOTE        0.5350  0.7392   0.6855 [0.6084 - 0.7617]             0.52
             Test_Specificity Test_Youden_J Test_F1
Sensitivity2           0.7912        0.4112  0.2183
Sensitivity5           0.8035        0.3835  0.2148
Sensitivity6           0.7109        0.4309  0.1962
Sensitivity4           0.6626        0.4226  0.1827
Sensitivity            0.6584        0.3384  0.1635
Sensitivity1           0.6348        0.3148  0.1549
Sensitivity3           0.7747        0.2747  0.1701
Sensitivity7           0.7623        0.2823  0.1694

09b Figure 4 — DAG of the primary model, with M3-vs-M6 structural check

library(bnlearn)
library(igraph)
library(ggraph)
library(ggplot2)
library(dplyr)

# ==============================================================================
# 1. Primary Model Setup & Structural Check
# ==============================================================================
dag_primary  <- dags[[PRIMARY]]
astr_primary <- arcs_l[[PRIMARY]]

cat("Plotting DAG for PRIMARY =", PRIMARY, "\n")

Plotting DAG for PRIMARY = M3: Continuous Nodes 

cat("Parents of stroke:", paste(parents(dag_primary, "stroke"), collapse = ", "), "\n")

Parents of stroke: hypertension, heart_disease, ever_married 

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

Children of stroke: age, avg_glucose_level 

# Structural check: Compare PRIMARY against M6
M6_NAME <- "M6: SICE + Continuous"
if (PRIMARY != M6_NAME && M6_NAME %in% names(dags)) {
  a_primary <- apply(arcs(dags[[PRIMARY]]), 1, paste, collapse = "->")
  a_m6      <- apply(arcs(dags[[M6_NAME]]), 1, paste, collapse = "->")
  only_primary <- setdiff(a_primary, a_m6)
  only_m6      <- setdiff(a_m6, a_primary)
  identical_structures <- length(only_primary) == 0 && length(only_m6) == 0

  cat("\n[Figure 4 check] Is", PRIMARY, "structurally identical to", M6_NAME, "?\n")
  cat("  Identical:", identical_structures, "\n")
  if (!identical_structures) {
    cat("  Arcs in", PRIMARY, "but not", M6_NAME, ":\n"); print(only_primary)
    cat("  Arcs in", M6_NAME, "but not", PRIMARY, ":\n"); print(only_m6)
  }
  write.csv(data.frame(only_in_primary = paste(only_primary, collapse = "; "),
                        only_in_m6      = paste(only_m6, collapse = "; "),
                        identical       = identical_structures),
            "v3_dag_M3_vs_M6_check.csv", row.names = FALSE)
}

[Figure 4 check] Is M3: Continuous Nodes structurally identical to M6: SICE + Continuous ?
  Identical: TRUE 

# ==============================================================================
# 2. Publication-Quality DAG Rendering via ggraph
# ==============================================================================
# Convert bnlearn DAG to igraph object
ig <- as.igraph(dag_primary)

# Map bootstrap arc strengths onto igraph edges
if (!is.null(astr_primary)) {
  str_df <- as.data.frame(astr_primary)
  e_mat  <- as_edgelist(ig)
  edge_strengths <- numeric(nrow(e_mat))
  
  for (i in seq_len(nrow(e_mat))) {
    u <- e_mat[i, 1]
    v <- e_mat[i, 2]
    match_idx <- which(str_df$from == u & str_df$to == v)
    if (length(match_idx) > 0) {
      edge_strengths[i] <- str_df$strength[match_idx[1]]
    } else {
      edge_strengths[i] <- 0.5
    }
  }
  E(ig)$strength <- edge_strengths
} else {
  E(ig)$strength <- 1.0
}

# Clean variable display labels
clean_labels <- c(
  "stroke"            = "Stroke\n(Target)",
  "age"               = "Age",
  "bmi"               = "BMI",
  "avg_glucose_level" = "Avg. Glucose",
  "hypertension"      = "Hypertension",
  "heart_disease"     = "Heart Disease",
  "ever_married"      = "Ever Married",
  "work_type"         = "Work Type",
  "smoking_status"    = "Smoking Status",
  "gender"            = "Gender",
  "Residence_type"    = "Residence Type"
)

# Define visual node roles
node_types <- c(
  "stroke"            = "Target Outcome",
  "age"               = "Continuous Predictor",
  "bmi"               = "Continuous Predictor",
  "avg_glucose_level" = "Continuous Predictor",
  "hypertension"      = "Categorical Predictor",
  "heart_disease"     = "Categorical Predictor",
  "ever_married"      = "Categorical Predictor",
  "work_type"         = "Categorical Predictor",
  "smoking_status"    = "Categorical Predictor",
  "gender"            = "Categorical Predictor",
  "Residence_type"    = "Categorical Predictor"
)

v_names <- V(ig)$name
V(ig)$label <- ifelse(v_names %in% names(clean_labels), clean_labels[v_names], v_names)
V(ig)$type  <- ifelse(v_names %in% names(node_types), node_types[v_names], "Categorical Predictor")
V(ig)$type  <- factor(V(ig)$type, levels = c("Target Outcome", "Continuous Predictor", "Categorical Predictor"))

# Plot using hierarchical Sugiyama layout
set.seed(if (exists("GLOBAL_SEED")) GLOBAL_SEED else 42)

p_dag <- ggraph(ig, layout = "sugiyama") +
  # Directed edges with strength-scaled thickness
  geom_edge_link(
    aes(edge_width = strength),
    arrow = arrow(length = unit(3.5, "mm"), type = "closed"),
    end_cap = circle(10, "mm"),
    start_cap = circle(10, "mm"),
    color = "#4A5568",
    alpha = 0.75
  ) +
  scale_edge_width_continuous(
    range = c(0.6, 2.2),
    name = "Bootstrap Arc Strength"
  ) +
  # Node styling
  geom_node_point(aes(color = type), size = 20) +
  geom_node_text(aes(label = label), fontface = "bold", size = 2.8, color = "#1A202C") +
  scale_color_manual(
    values = c(
      "Target Outcome"        = "#E53E3E", # Red/Coral fill
      "Continuous Predictor"  = "#3182CE", # Soft Blue fill
      "Categorical Predictor" = "#E2E8F0"  # Light Slate Gray fill
    ),
    name = "Node Category"
  ) +
  theme_void() +
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.title = element_text(face = "bold", size = 10),
    legend.text = element_text(size = 9),
    plot.margin = margin(15, 15, 15, 15)
  )

# Export clean 300 DPI image
ggsave("v3_DAG_primary.png", plot = p_dag, width = 10, height = 7.5, dpi = 300)
cat("Saved publication-ready v3_DAG_primary.png\n")

Saved publication-ready v3_DAG_primary.png

# ==============================================================================
# 3. Model Diagnostics & Logs
# ==============================================================================
cat("\nModel string:\n", modelstring(dag_primary), "\n")

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("Full arc list (from -> to):\n"); print(arcs(dag_primary))

Full 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,] "bmi"           "avg_glucose_level"
[14,] "heart_disease" "ever_married"     
[15,] "age"           "bmi"              
[16,] "gender"        "heart_disease"    
[17,] "ever_married"  "stroke"           

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

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

10 ROC figure

cb <- c("#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7")
rp <- roc_list; names(rp) <- gsub(" \\(.*\\)","",names(rp))
g <- pROC::ggroc(rp, legacy.axes=TRUE, size=1) +
  geom_abline(intercept=0, slope=1, colour="grey50", linetype="dashed") +
  scale_color_manual(values=cb) +
  labs(x="False Positive Rate (1 - Specificity)", y="True Positive Rate (Sensitivity)") +
  theme_minimal() + theme(legend.title=element_blank())
ggsave("v3_ROC.png", g, width=8, height=6, dpi=300); print(g)

11 Calibration — reference model, weak calibration, prior correction

#A prevalence-only model attains Brier ~= p(1-p) ~= 0.047 here, so raw
# Brier alone cannot demonstrate calibration. BSS, intercept and slope are added.
calib_metrics <- function(y,p,eps=1e-6){
  y <- as.numeric(as.character(y)); p <- pmin(pmax(p,eps),1-eps)
  brier <- mean((p-y)^2); ref <- mean((mean(y)-y)^2); lp <- qlogis(p)
  data.frame(Brier=round(brier,4), Brier_ref=round(ref,4),
             BSS=round(1-brier/ref,4),
             Calib_intercept=round(unname(coef(glm(y~offset(lp),family=binomial))[1]),3),
             Calib_slope=round(unname(coef(glm(y~lp,family=binomial))[2]),3))
}
prior_correct <- function(p, pi_train=0.5, pi_pop=PI_POP, eps=1e-6){
  p <- pmin(pmax(p,eps),1-eps)
  plogis(qlogis(p) - log((pi_train/(1-pi_train))/(pi_pop/(1-pi_pop))))
}
SMOTE_ARMS <- c("M4","M5","M7","M8")
calib <- map_dfr(names(roc_list), function(m){
  y <- roc_list[[m]]$response; p <- roc_list[[m]]$predictor
  out <- cbind(Model=m, Correction="none", calib_metrics(y,p))
  if (any(sapply(SMOTE_ARMS, function(a) startsWith(m,a))))
    out <- bind_rows(out, cbind(Model=m, Correction="prior-corrected",
                                calib_metrics(y, prior_correct(p))))
  out })
print(calib); write.csv(calib,"v3_calibration.csv",row.names=FALSE)

                           Model      Correction  Brier Brier_ref     BSS Calib_intercept Calib_slope
1                   M1: Baseline            none 0.0504    0.0465 -0.0824          -0.206       0.257
2          M2: SICE + Discretize            none 0.0509    0.0465 -0.0941          -0.217       0.240
3           M3: Continuous Nodes            none 0.0561    0.0465 -0.2049          -0.643       0.344
4          M4: SMOTE on Baseline            none 0.1735    0.0465 -2.7279          -3.012       0.153
5          M4: SMOTE on Baseline prior-corrected 0.0489    0.0465 -0.0519          -0.041       0.174
6  M5: ALL (SICE + Cont + SMOTE)            none 0.1756    0.0465 -2.7739          -3.282       0.174
7  M5: ALL (SICE + Cont + SMOTE) prior-corrected 0.0542    0.0465 -0.1651          -0.311       0.193
8          M6: SICE + Continuous            none 0.0564    0.0465 -0.2132          -0.651       0.345
9         M7: Continuous + SMOTE            none 0.1768    0.0465 -2.8007          -3.401       0.165
10        M7: Continuous + SMOTE prior-corrected 0.0552    0.0465 -0.1859          -0.430       0.188
11              M8: SICE + SMOTE            none 0.1750    0.0465 -2.7611          -3.049       0.139
12              M8: SICE + SMOTE prior-corrected 0.0493    0.0465 -0.0593          -0.077       0.169

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))
  gg <- if (length(br)<3) cut(p, breaks=unique(c(-Inf,sort(unique(p)))), include.lowest=TRUE)
        else cut(p, breaks=br, include.lowest=TRUE)
  na.omit(data.frame(model=m, pred=tapply(p,gg,mean), obs=tapply(y,gg,mean)))
}
focus <- intersect(c("M1: Baseline","M3: Continuous Nodes","M6: SICE + Continuous",
                     "M5: ALL (SICE + Cont + SMOTE)"), names(roc_list))
cp <- ggplot(bind_rows(lapply(focus,calib_curve)), aes(pred,obs,colour=model)) +
  geom_abline(slope=1,intercept=0,linetype="dashed",colour="grey50") +
  geom_hline(yintercept=PI_POP, linetype="dotted", colour="red") +
  geom_line()+geom_point() +
  labs(x="Mean predicted probability", y="Observed stroke fraction",
       caption="Dotted red line: prevalence-only reference") + theme_minimal()
ggsave("v3_calibration.png", cp, width=7, height=6, dpi=300); print(cp)

12 Comparison with alternative models (transparent + opaque) and SHAP

bt <- train_m3; bv <- val_m3; be <- test_m3

# 1. Logistic Regression
lr <- glm(stroke ~ ., data = bt, family = binomial)

# FIX: Save the actual ROC object to 'roc_lr'
roc_lr <- roc(be$stroke, predict(lr, be, type = "response"), levels = c("0", "1"), quiet = TRUE)
lr_auc <- as.numeric(auc(roc_lr))

# 2. XGBoost
mm <- function(d) model.matrix(stroke ~ . -1, data = d)
ft <- Reduce(intersect, list(colnames(mm(bt)), colnames(mm(bv)), colnames(mm(be))))

dtr <- xgb.DMatrix(mm(bt)[, ft, drop = FALSE], label = as.numeric(as.character(bt$stroke)))
dva <- xgb.DMatrix(mm(bv)[, ft, drop = FALSE], label = as.numeric(as.character(bv$stroke)))
dte <- xgb.DMatrix(mm(be)[, ft, drop = FALSE], label = as.numeric(as.character(be$stroke)))

SPW <- sum(bt$stroke == "0") / sum(bt$stroke == "1")

XGB_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 = SPW
)

set.seed(GLOBAL_SEED)
xgb <- xgb.train(
  params = XGB_PARAMS, 
  data = dtr, 
  nrounds = 1000,
  watchlist = list(val = dva), 
  early_stopping_rounds = 25, 
  verbose = 0
)

# FIX: Save the actual ROC object to 'roc_xgb'
roc_xgb <- roc(be$stroke, predict(xgb, dte), levels = c("0", "1"), quiet = TRUE)
xgb_auc <- as.numeric(auc(roc_xgb))

# Hyperparameter output for Methods
cat("\n===== XGBoost specification (paste into Methods) =====\n")

===== XGBoost specification (paste into Methods) =====

cat(sprintf("objective=%s; eval_metric=%s; max_depth=%d; learning_rate=%.2f;\n",
    XGB_PARAMS$objective, XGB_PARAMS$eval_metric, XGB_PARAMS$max_depth, XGB_PARAMS$learning_rate))

objective=binary:logistic; eval_metric=auc; max_depth=3; learning_rate=0.05;

cat(sprintf("subsample=%.1f; colsample_bytree=%.1f; scale_pos_weight=%.3f;\n",
    XGB_PARAMS$subsample, XGB_PARAMS$colsample_bytree, SPW))

subsample=0.8; colsample_bytree=0.8; scale_pos_weight=19.570;

cat(sprintf("early_stopping_rounds=25; n_rounds_selected=%d; seed=%d\n",
    xgb$best_iteration, GLOBAL_SEED))

bench <- data.frame(
  Model = c(paste0("Bayesian network (", PRIMARY, ")"),
            "Logistic regression (transparent)",
            "XGBoost (opaque)"),
  Test_AUC = round(c(leaderboard$Test_AUC[leaderboard$Model == PRIMARY], lr_auc, xgb_auc), 4)
)
print(bench)

                                    Model Test_AUC
1 Bayesian network (M3: Continuous Nodes)   0.8109
2       Logistic regression (transparent)   0.8188
3                        XGBoost (opaque)   0.8222

write.csv(bench, "v3_benchmark.csv", row.names = FALSE)

# [R3.15] SHAP for the opaque benchmark
sv  <- predict(xgb, dte, predcontrib = TRUE)
imp <- sort(colMeans(abs(sv[, setdiff(colnames(sv), "BIAS"), drop = FALSE])), decreasing = TRUE)
print(round(head(imp, 10), 5))

                       age                        bmi          avg_glucose_level               genderFemale 
                   0.85756                    0.15479                    0.10572                    0.02464 
             hypertension1          work_typeGovt_job       smoking_statussmokes                (Intercept) 
                   0.02349                    0.01622                    0.01375                    0.01166 
smoking_statusnever smoked     work_typeSelf-employed 
                   0.01138                    0.01059 

write.csv(data.frame(feature = names(imp), mean_abs_shap = round(imp, 5)),
          "v3_shap.csv", row.names = FALSE)

13 Hypothesis testing — DeLong and McNemar (complete set)

thr <- function(m) leaderboard$Val_Threshold[leaderboard$Model==m]

mcn <- 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); i <- which(y==1)
  p <- suppressWarnings(mcnemar.test(table(factor(pa[i],0:1),factor(pb[i],0:1)))$p.value)
  cat(sprintf("McNemar  %-30s vs %-30s p=%.4f\n",A,B,p)); p 
}

dl <- function(A,B){
  t <- roc.test(roc_list[[A]], roc_list[[B]], method="delong")
  cat(sprintf("DeLong   %-30s vs %-30s Z=%6.3f p=%.4f\n",A,B,
      as.numeric(t$statistic),t$p.value)); c(Z=as.numeric(t$statistic),p=t$p.value) 
}

M5 <- "M5: ALL (SICE + Cont + SMOTE)"

cat("\n== H1 topology ==\n"); dl("M3: Continuous Nodes","M1: Baseline"); dl("M6: SICE + Continuous","M2: SICE + Discretize")

== H1 topology ==
DeLong   M3: Continuous Nodes           vs M1: Baseline                   Z= 2.402 p=0.0163
         Z          p 
2.40204355 0.01630377 
DeLong   M6: SICE + Continuous          vs M2: SICE + Discretize          Z= 2.630 p=0.0085
          Z           p 
2.629591796 0.008548745 

mcn("M3: Continuous Nodes","M1: Baseline"); mcn("M6: SICE + Continuous","M2: SICE + Discretize")

McNemar  M3: Continuous Nodes           vs M1: Baseline                   p=0.6276
[1] 0.6276258
McNemar  M6: SICE + Continuous          vs M2: SICE + Discretize          p=0.3827
[1] 0.3827331

cat("\n== H2 imputation ==\n"); dl("M2: SICE + Discretize","M1: Baseline"); dl("M6: SICE + Continuous","M3: Continuous Nodes")

== H2 imputation ==
DeLong   M2: SICE + Discretize          vs M1: Baseline                   Z=-1.399 p=0.1618
         Z          p 
-1.3989748  0.1618205 
DeLong   M6: SICE + Continuous          vs M3: Continuous Nodes           Z=-0.297 p=0.7667
         Z          p 
-0.2967460  0.7666604 

cat("\n== H3 balancing ==\n"); dl("M4: SMOTE on Baseline","M1: Baseline"); dl(M5,"M6: SICE + Continuous")

== H3 balancing ==
DeLong   M4: SMOTE on Baseline          vs M1: Baseline                   Z=-1.893 p=0.0584
          Z           p 
-1.89266368  0.05840261 
DeLong   M5: ALL (SICE + Cont + SMOTE)  vs M6: SICE + Continuous          Z=-2.060 p=0.0394
          Z           p 
-2.05954502  0.03944206 

mcn("M4: SMOTE on Baseline","M1: Baseline"); mcn("M8: SICE + SMOTE","M2: SICE + Discretize"); mcn(M5,"M6: SICE + Continuous")

McNemar  M4: SMOTE on Baseline          vs M1: Baseline                   p=0.0077
[1] 0.007660761
McNemar  M8: SICE + SMOTE               vs M2: SICE + Discretize          p=0.0133
[1] 0.01332833
McNemar  M5: ALL (SICE + Cont + SMOTE)  vs M6: SICE + Continuous          p=0.0809
[1] 0.0808556

cat("\n== H4 combinations ==\n"); dl(M5,"M1: Baseline"); dl("M6: SICE + Continuous","M1: Baseline")

== H4 combinations ==
DeLong   M5: ALL (SICE + Cont + SMOTE)  vs M1: Baseline                   Z= 0.338 p=0.7352
        Z         p 
0.3382692 0.7351604 
DeLong   M6: SICE + Continuous          vs M1: Baseline                   Z= 2.338 p=0.0194
         Z          p 
2.33845878 0.01936346 

# ==============================================================================
# Benchmark Model Comparisons (BN vs Logistic Regression & XGBoost)
# ==============================================================================
cat("\n== Benchmark Comparisons ==\n")

== Benchmark Comparisons ==

# roc_bn MUST be the validation-selected PRIMARY model (set in Segment 09), not M5.
# M5 is the fully-combined SICE+continuous+SMOTE arm and is a much weaker
# configuration (test AUC 0.7381) than the model actually reported in the
# manuscript's leaderboard (PRIMARY = M3, test AUC 0.8109). Comparing benchmarks
# against M5 silently answers a different question -- it does not error, it just
# produces the wrong p-values -- which is what happened in the previous run.
roc_bn <- roc_list[[PRIMARY]]

t_lr  <- roc.test(roc_bn, roc_lr,  method = "delong")
t_xgb <- roc.test(roc_bn, roc_xgb, method = "delong")

cat(sprintf("DeLong   %-30s vs %-30s Z=%6.3f p=%.4f\n",
            PRIMARY, "Logistic regression", as.numeric(t_lr$statistic),  t_lr$p.value))

DeLong   M3: Continuous Nodes           vs Logistic regression            Z=-0.735 p=0.4622

cat(sprintf("DeLong   %-30s vs %-30s Z=%6.3f p=%.4f\n",
            PRIMARY, "XGBoost",             as.numeric(t_xgb$statistic), t_xgb$p.value))

DeLong   M3: Continuous Nodes           vs XGBoost                        Z=-1.010 p=0.3126

14 Robustness — all arms, repeated stratified resampling

eval_light <- function(tr, te, disc){
  # 1. Dynamically synchronize factor levels between training (tr) and test (te) sets
  for (col in names(tr)) {
    if (is.factor(tr[[col]])) {
      # Remove unobserved levels from tr (fixes check.data warning)
      tr[[col]] <- droplevels(tr[[col]]) 
      # Force te to match tr's exact factor levels (fixes level mismatch skip error)
      te[[col]] <- factor(te[[col]], levels = levels(tr[[col]])) 
    }
  }
  
  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)))
}

safe <- function(e) tryCatch(e, error = function(x){ cat("  skip:", x$message, "\n"); NA_real_ })

rob <- data.frame()
for (s in 1:N_REPEATS) {
  q <- split_stroke(stroke_data, seed = 1000 + s)
  q <- harmonise_levels(q$train, q$val, q$test)
  
  m <- median(as.numeric(as.character(q$train$bmi)), na.rm = TRUE)
  ce <- impute_sice(q$train, q$val, q$test, seed = 1000 + s)
  
  d1 <- impute_baseline(q$train, m); d1t <- impute_baseline(q$test, m)
  d2 <- impute_baseline(ce$train, m); d2t <- impute_baseline(ce$test, m)
  c3 <- prep_continuous(q$train, m); c3t <- prep_continuous(q$test, m)
  c6 <- prep_continuous(ce$train, m); c6t <- prep_continuous(ce$test, m)
  
  rob <- rbind(rob, data.frame(seed = s,
    M1 = safe(eval_light(d1, d1t, TRUE)),
    M2 = safe(eval_light(d2, d2t, TRUE)),
    M3 = safe(eval_light(c3, c3t, FALSE)),
    M4 = safe(eval_light(apply_smotenc(d1), d1t, TRUE)),
    M5 = safe(eval_light(apply_smotenc(c6), c6t, FALSE)),
    M6 = safe(eval_light(c6, c6t, FALSE)),
    M7 = safe(eval_light(apply_smotenc(c3), c3t, FALSE)),
    M8 = safe(eval_light(apply_smotenc(d2), d2t, TRUE))
  ))
  cat(sprintf("seed %d/%d\n", s, N_REPEATS))
}

seed 1/30
seed 2/30
seed 3/30
seed 4/30
seed 5/30
seed 6/30
seed 7/30
seed 8/30
seed 9/30
seed 10/30
seed 11/30
seed 12/30
seed 13/30
seed 14/30
seed 15/30
seed 16/30
seed 17/30
seed 18/30
seed 19/30
seed 20/30
seed 21/30
seed 22/30
seed 23/30
seed 24/30
seed 25/30
seed 26/30
seed 27/30
seed 28/30
seed 29/30
seed 30/30

cat("\n ANALYSABLE RESAMPLES PER ARM:\n")

 ANALYSABLE RESAMPLES PER ARM:

print(colSums(!is.na(rob[-1])))

M1 M2 M3 M4 M5 M6 M7 M8 
30 30 30 30 30 30 30 30 

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

rob_summary <- rob %>% 
  pivot_longer(-seed, names_to = "arm", values_to = "auc") %>%
  group_by(arm) %>% 
  summarise(
    n_analysable = sum(!is.na(auc)),
    AUC_mean = round(mean(auc, na.rm = TRUE), 4),
    AUC_sd = round(sd(auc,na.rm = TRUE), 4),
    .groups = "drop"
  ) %>% 
  arrange(desc(AUC_mean))

print(rob_summary)

# A tibble: 8 × 4
  arm   n_analysable AUC_mean AUC_sd
  <chr>        <int>    <dbl>  <dbl>
1 M3              30    0.820 0.0304
2 M6              30    0.817 0.0297
3 M1              30    0.772 0.037 
4 M2              30    0.772 0.0392
5 M4              30    0.758 0.0371
6 M8              30    0.756 0.0306
7 M5              30    0.746 0.0406
8 M7              30    0.746 0.0409

write.csv(rob_summary, "v3_robustness_summary.csv", row.names = FALSE)

15 Factorial analysis of the 2×2×2 design

# Main effects AND interactions. If no interaction is significant, the
# word "synergy" must be removed from Hypothesis 4 and the Conclusion.
d <- rob %>% pivot_longer(-seed,names_to="arm",values_to="auc") %>% filter(!is.na(auc)) %>%
  mutate(topology  =factor(if_else(arm %in% c("M3","M5","M6","M7"),"continuous","discrete"),
                           levels=c("discrete","continuous")),
         imputation=factor(if_else(arm %in% c("M2","M5","M6","M8"),"sice","median"),
                           levels=c("median","sice")),
         balancing =factor(if_else(arm %in% c("M4","M5","M7","M8"),"smote","none"),
                           levels=c("none","smote")), seed=factor(seed))
fm <- lmer(auc ~ topology*imputation*balancing + (1|seed), data=d, REML=TRUE)
cat("\n--- Type III fixed effects [R1.3a] ---\n"); print(anova(fm, type=3))

--- Type III fixed effects [R1.3a] ---
Type III Analysis of Variance Table with Satterthwaite's method
                                Sum Sq  Mean Sq NumDF DenDF  F value    Pr(>F)    
topology                      0.018368 0.018368     1   203  29.1913 1.827e-07 ***
imputation                    0.000068 0.000068     1   203   0.1078    0.7430    
balancing                     0.112374 0.112374     1   203 178.5938 < 2.2e-16 ***
topology:imputation           0.000000 0.000000     1   203   0.0000    0.9946    
topology:balancing            0.050988 0.050988     1   203  81.0342 < 2.2e-16 ***
imputation:balancing          0.000005 0.000005     1   203   0.0077    0.9303    
topology:imputation:balancing 0.000097 0.000097     1   203   0.1539    0.6952    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

cat("\n--- Estimates with 95% CI ---\n")

--- Estimates with 95% CI ---

print(round(cbind(Estimate=fixef(fm), confint(fm,method="Wald")[names(fixef(fm)),]),4))

                                                 Estimate   2.5 %  97.5 %
(Intercept)                                        0.7716  0.7588  0.7845
topologycontinuous                                 0.0479  0.0352  0.0606
imputationsice                                    -0.0001 -0.0127  0.0126
balancingsmote                                    -0.0131 -0.0258 -0.0004
topologycontinuous:imputationsice                 -0.0026 -0.0205  0.0154
topologycontinuous:balancingsmote                 -0.0608 -0.0788 -0.0429
imputationsice:balancingsmote                     -0.0020 -0.0199  0.0160
topologycontinuous:imputationsice:balancingsmote   0.0051 -0.0203  0.0305

for (v in c("topology","imputation","balancing"))
  print(d %>% group_by(.data[[v]]) %>% summarise(mean_auc=round(mean(auc),4), sd=round(sd(auc),4)))

# A tibble: 2 × 3
  topology   mean_auc     sd
  <fct>         <dbl>  <dbl>
1 discrete      0.765 0.0364
2 continuous    0.782 0.0507
# A tibble: 2 × 3
  imputation mean_auc     sd
  <fct>         <dbl>  <dbl>
1 median        0.774 0.0457
2 sice          0.773 0.0443
# A tibble: 2 × 3
  balancing mean_auc     sd
  <fct>        <dbl>  <dbl>
1 none         0.795 0.0412
2 smote        0.752 0.0375

w <- rob %>% filter(!is.na(M3),!is.na(M1))
cat("\nPaired confirmatory (M3 vs M1):\n"); print(t.test(w$M3,w$M1,paired=TRUE)); print(wilcox.test(w$M3,w$M1,paired=TRUE))

Paired confirmatory (M3 vs M1):

    Paired t-test

data:  w$M3 and w$M1
t = 6.8021, df = 29, p-value = 1.812e-07
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.03352577 0.06235489
sample estimates:
mean difference 
     0.04794033 


    Wilcoxon signed rank exact test

data:  w$M3 and w$M1
V = 441, p-value = 1.419e-06
alternative hypothesis: true location shift is not equal to 0

16 Sparse CPT reliability and consensus network

# ---- Fit the bootstrap-averaged CONSENSUS network and evaluate it ----
pr <- runs[[PRIMARY]]
astr <- pr$arc_strength
cons <- averaged.network(astr)
cons_directed <- TRUE
cons <- tryCatch(cextend(cons), error = function(e) {
  cons_directed <<- FALSE
  averaged.network(astr, threshold = 0.85)
})

cat("\n[R1.5a] Consensus network arcs:\n"); print(arcs(cons))

[R1.5a] Consensus network arcs:
      from            to                 
 [1,] "age"           "avg_glucose_level"
 [2,] "age"           "bmi"              
 [3,] "hypertension"  "stroke"           
 [4,] "heart_disease" "stroke"           
 [5,] "ever_married"  "age"              
 [6,] "ever_married"  "heart_disease"    
 [7,] "ever_married"  "avg_glucose_level"
 [8,] "ever_married"  "stroke"           
 [9,] "work_type"     "age"              
[10,] "work_type"     "hypertension"     
[11,] "work_type"     "ever_married"     
[12,] "work_type"     "bmi"              
[13,] "work_type"     "smoking_status"   
[14,] "bmi"           "avg_glucose_level"
[15,] "stroke"        "age"              
[16,] "stroke"        "avg_glucose_level"

cat("\nArcs in single HC fit but not consensus:\n")

Arcs in single HC fit but not consensus:

h <- apply(arcs(pr$structure),1,paste,collapse="->"); c2 <- apply(arcs(cons),1,paste,collapse="->")
print(setdiff(h,c2)); cat("Arcs in consensus but not HC fit:\n"); print(setdiff(c2,h))

[1] "ever_married->work_type"     "heart_disease->ever_married" "gender->heart_disease"      
Arcs in consensus but not HC fit:
[1] "ever_married->heart_disease" "work_type->ever_married"    

cons_auc <- tryCatch({
  cf <- if (pr$is_disc) bn.fit(cons, pr$train, method="bayes", iss=10) else bn.fit(cons, pr$train)
  gp <- if (pr$is_disc) {
    attr(predict(cf,node="stroke",data=pr$test,method="bayes-lw",prob=TRUE),"prob")["1",]
  } else {
    ec <- setdiff(names(pr$test),"stroke")
    vapply(seq_len(nrow(pr$test)),function(i){
      p<-suppressWarnings(cpquery(cf,event=(stroke=="1"),evidence=as.list(pr$test[i,ec]),
         method="lw",n=N_PARTICLES_PATIENT)); if(is.na(p)) 0 else p},numeric(1)) }
  round(as.numeric(auc(roc(pr$test$stroke, gp, levels=c("0","1"), quiet=TRUE))),4)
}, error=function(e){cat("consensus fit failed:",e$message,"\n"); NA})
cat(sprintf("\n[R1.5a] Consensus-network test AUC = %s (single HC fit = %.4f)\n",
            cons_auc, leaderboard$Test_AUC[leaderboard$Model==PRIMARY]))

[R1.5a] Consensus-network test AUC = 0.8115 (single HC fit = 0.8109)

write.csv(astr %>% filter(strength>0.5, direction>=0.5) %>% arrange(desc(strength)),
          "v3_arc_strengths.csv", row.names=FALSE)

# ---- CPT counts + Jeffreys intervals + shrinkage ----
cpt_tbl <- train_m6 %>% count(hypertension, heart_disease, ever_married, stroke) %>%
  pivot_wider(names_from=stroke, values_from=n, values_fill=0) %>%
  rename(n0=`0`, n1=`1`) %>%
  mutate(n=n0+n1, p_hat=round(100*n1/n,2),
         CrI_low =round(100*qbeta(.025,n1+.5,n0+.5),2),
         CrI_high=round(100*qbeta(.975,n1+.5,n0+.5),2)) %>% arrange(desc(p_hat))
print(cpt_tbl); write.csv(cpt_tbl,"v3_cpt_intervals.csv",row.names=FALSE)

# A tibble: 8 × 9
  hypertension heart_disease ever_married    n0    n1     n p_hat CrI_low CrI_high
  <fct>        <fct>         <fct>        <int> <int> <int> <dbl>   <dbl>    <dbl>
1 1            1             No               1     3     4 75      28.4     97.2 
2 0            1             No              11     2    13 15.4     3.34    40.9 
3 1            0             No              26     4    30 13.3     4.67    28.6 
4 1            1             Yes             26     4    30 13.3     4.67    28.6 
5 0            1             Yes             94    14   108 13.0     7.63    20.2 
6 1            0             Yes            213    26   239 10.9     7.4     15.3 
7 0            0             Yes           1551    84  1635  5.14    4.15     6.29
8 0            0             No             994    12  1006  1.19    0.65     2.01

p_marg <- mean(train_m6$stroke=="1")
iss_tbl <- map_dfr(c(1,5,10,25,50,100), function(k)
  cpt_tbl %>% transmute(hypertension,heart_disease,ever_married,n,iss=k,
                        p_smoothed=round(100*(n1+k*p_marg)/(n+k),2)))
write.csv(iss_tbl,"v3_cpt_iss_sensitivity.csv",row.names=FALSE)
print(iss_tbl %>% filter(n<=10))

# A tibble: 6 × 6
  hypertension heart_disease ever_married     n   iss p_smoothed
  <fct>        <fct>         <fct>        <int> <dbl>      <dbl>
1 1            1             No               4     1      61.0 
2 1            1             No               4     5      36.0 
3 1            1             No               4    10      24.9 
4 1            1             No               4    25      14.5 
5 1            1             No               4    50      10.1 
6 1            1             No               4   100       7.56

# ---- [R3.8] Does the DISCRETE arm recover age -> stroke? ----
a1 <- arcs(dags[["M1: Baseline"]])
cat("\n[R3.8] Discrete-arm arcs involving age:\n"); print(a1[apply(a1,1,function(r) any(r=="age")),,drop=FALSE])

[R3.8] Discrete-arm arcs involving age:
     from  to             
[1,] "age" "stroke"       
[2,] "age" "work_type"    
[3,] "age" "ever_married" 
[4,] "age" "hypertension" 
[5,] "age" "heart_disease"

17 Conditional-risk scenarios (NOT counterfactual interventions)

# [R1.1] Observational conditioning, not do(). No treatment-effect language.
net <- fits[[PRIMARY]]; set.seed(2026)
pA <- cpquery(net, event=(stroke=="1"),
              evidence=(age<=50 & hypertension=="0" & heart_disease=="0"), n=N_PARTICLES_QUERY)
pB <- cpquery(net, event=(stroke=="1"),
              evidence=(age>=65 & hypertension=="1" & heart_disease=="1" & avg_glucose_level>=150), n=N_PARTICLES_QUERY)
pC <- cpquery(net, event=(stroke=="1"),
              evidence=(age>=65 & hypertension=="1" & heart_disease=="1" & avg_glucose_level<=100), n=N_PARTICLES_QUERY)
scen <- data.frame(Profile=c("A: age<=50, no HTN, no HD",
                             "B: age>=65, HTN, HD, glucose>=150",
                             "C: age>=65, HTN, HD, glucose<=100"),
                   Conditional_probability=sprintf("%.2f%%",100*c(pA,pB,pC)))
print(scen); write.csv(scen,"v3_conditional_risk.csv",row.names=FALSE)

                            Profile Conditional_probability
1         A: age<=50, no HTN, no HD                   0.74%
2 B: age>=65, HTN, HD, glucose>=150                  46.13%
3 C: age>=65, HTN, HD, glucose<=100                  29.89%

18 Session information

writeLines(capture.output(sessionInfo()), "v2_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] ggraph_2.2.2        igraph_2.2.1        vcd_1.4-13          xgboost_3.2.1.1     lmerTest_3.2-1     
 [6] lme4_1.1-38         Matrix_1.7-4        themis_1.0.3        recipes_1.3.1       mice_3.19.0        
[11] pROC_1.19.0.1       caret_7.0-1         lattice_0.22-7      bnlearn_5.1         lubridate_1.9.4    
[16] forcats_1.0.1       stringr_1.6.0       dplyr_1.1.4         purrr_1.2.0         readr_2.1.6        
[21] tidyr_1.3.1         tibble_3.3.0        ggplot2_4.0.1       tidyverse_2.0.0     graph_1.88.1       
[26] BiocGenerics_0.56.0 generics_0.1.4     

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

19 RESULTS PACK — every number the Results section needs

# Run last. Produces a single console block ordered to match the manuscript's
# Results subsections 3.1 - 3.9, so numbers can be transcribed directly.
line <- function(x) cat(strrep("=",72),"\n",x,"\n",strrep("=",72),"\n",sep="")

line("3.0  COHORT AND PARTITIONS")

========================================================================
3.0  COHORT AND PARTITIONS
========================================================================

cat(sprintf("N = %d; prevalence = %.2f%%; missing BMI = %.1f%%\n",
            N_TOTAL, 100*PI_POP, 100*mean(is.na(stroke_data$bmi))))

N = 5109; prevalence = 4.87%; missing BMI = 3.9%

print(split_tbl)

# A tibble: 3 × 4
  Partition      n Events Prevalence
  <chr>      <int>  <int> <chr>     
1 Training    3065    149 4.86%     
2 Validation  1022     50 4.89%     
3 Test        1022     50 4.89%     

line("3.1  TRAINING-ONLY ASSOCIATION SCREEN (Supplementary)")

========================================================================
3.1  TRAINING-ONLY ASSOCIATION SCREEN (Supplementary)
========================================================================

print(assoc_tbl)

            Variable        Type  Assoc
1                age  continuous 0.2410
2  avg_glucose_level  continuous 0.1318
3       hypertension categorical 0.1132
4      heart_disease categorical 0.1071
5       ever_married categorical 0.0964
6          work_type categorical 0.0926
7     smoking_status categorical 0.0750
8                bmi  continuous 0.0356
9     Residence_type categorical 0.0250
10            gender categorical 0.0008

line("3.2  STRUCTURE — primary model, consensus comparison, arc direction")

========================================================================
3.2  STRUCTURE — primary model, consensus comparison, arc direction
========================================================================

cat("PRIMARY MODEL (validation-selected):", PRIMARY, "\n")

PRIMARY MODEL (validation-selected): M3: Continuous Nodes 

cat("Parents of stroke:", paste(parents(fits[[PRIMARY]],"stroke"), collapse=", "), "\n")

Parents of stroke: hypertension, heart_disease, ever_married 

cat("Children of stroke:", paste(children(fits[[PRIMARY]],"stroke"), collapse=", "), "\n")

Children of stroke: age, avg_glucose_level 

cat(sprintf("Consensus-network test AUC = %s vs single HC = %.4f\n",
            cons_auc, leaderboard$Test_AUC[leaderboard$Model==PRIMARY]))

Consensus-network test AUC = 0.8115 vs single HC = 0.8109

line("3.3  CPT RELIABILITY")

========================================================================
3.3  CPT RELIABILITY
========================================================================

print(cpt_tbl); cat("\nShrinkage for sparse cells:\n"); print(iss_tbl %>% filter(n<=10))

# A tibble: 8 × 9
  hypertension heart_disease ever_married    n0    n1     n p_hat CrI_low CrI_high
  <fct>        <fct>         <fct>        <int> <int> <int> <dbl>   <dbl>    <dbl>
1 1            1             No               1     3     4 75      28.4     97.2 
2 0            1             No              11     2    13 15.4     3.34    40.9 
3 1            0             No              26     4    30 13.3     4.67    28.6 
4 1            1             Yes             26     4    30 13.3     4.67    28.6 
5 0            1             Yes             94    14   108 13.0     7.63    20.2 
6 1            0             Yes            213    26   239 10.9     7.4     15.3 
7 0            0             Yes           1551    84  1635  5.14    4.15     6.29
8 0            0             No             994    12  1006  1.19    0.65     2.01

Shrinkage for sparse cells:
# A tibble: 6 × 6
  hypertension heart_disease ever_married     n   iss p_smoothed
  <fct>        <fct>         <fct>        <int> <dbl>      <dbl>
1 1            1             No               4     1      61.0 
2 1            1             No               4     5      36.0 
3 1            1             No               4    10      24.9 
4 1            1             No               4    25      14.5 
5 1            1             No               4    50      10.1 
6 1            1             No               4   100       7.56

line("3.4  CONDITIONAL-RISK SCENARIOS")

========================================================================
3.4  CONDITIONAL-RISK SCENARIOS
========================================================================

print(scen)

                            Profile Conditional_probability
1         A: age<=50, no HTN, no HD                   0.74%
2 B: age>=65, HTN, HD, glucose>=150                  46.13%
3 C: age>=65, HTN, HD, glucose<=100                  29.89%

line("3.5  PERFORMANCE LEADERBOARD")

========================================================================
3.5  PERFORMANCE LEADERBOARD
========================================================================

print(leaderboard)

                                     Model Val_Threshold Val_AUC Test_AUC    Test_AUC_95_CI Test_Sensitivity
Sensitivity2          M3: Continuous Nodes        0.0857  0.8642   0.8109 [0.7612 - 0.8533]             0.62
Sensitivity5         M6: SICE + Continuous        0.0951  0.8683   0.8094 [0.7615 - 0.8522]             0.58
Sensitivity6        M7: Continuous + SMOTE        0.4303  0.7602   0.7456 [0.6672 - 0.8126]             0.72
Sensitivity4 M5: ALL (SICE + Cont + SMOTE)        0.3022  0.7500   0.7381 [0.6688 - 0.8048]             0.76
Sensitivity                   M1: Baseline        0.0350  0.7710   0.7235 [0.6542 - 0.7907]             0.68
Sensitivity1         M2: SICE + Discretize        0.0330  0.7774   0.7108 [0.6376 - 0.7805]             0.68
Sensitivity3         M4: SMOTE on Baseline        0.5530  0.7509   0.6903 [0.6128 - 0.7630]             0.50
Sensitivity7              M8: SICE + SMOTE        0.5350  0.7392   0.6855 [0.6084 - 0.7617]             0.52
             Test_Specificity Test_Youden_J Test_F1
Sensitivity2           0.7912        0.4112  0.2183
Sensitivity5           0.8035        0.3835  0.2148
Sensitivity6           0.7109        0.4309  0.1962
Sensitivity4           0.6626        0.4226  0.1827
Sensitivity            0.6584        0.3384  0.1635
Sensitivity1           0.6348        0.3148  0.1549
Sensitivity3           0.7747        0.2747  0.1701
Sensitivity7           0.7623        0.2823  0.1694

line("3.6  HYPOTHESIS TESTS")

========================================================================
3.6  HYPOTHESIS TESTS
========================================================================

cat("(see Segment 13 console output; transcribe Z and p for H1-H4)\n")

(see Segment 13 console output; transcribe Z and p for H1-H4)

line("3.7  ROBUSTNESS")

========================================================================
3.7  ROBUSTNESS
========================================================================

print(rob_summary)

# A tibble: 8 × 4
  arm   n_analysable AUC_mean AUC_sd
  <chr>        <int>    <dbl>  <dbl>
1 M3              30    0.820 0.0304
2 M6              30    0.817 0.0297
3 M1              30    0.772 0.037 
4 M2              30    0.772 0.0392
5 M4              30    0.758 0.0371
6 M8              30    0.756 0.0306
7 M5              30    0.746 0.0406
8 M7              30    0.746 0.0409

cat("\nAnalysable resamples per arm:\n"); print(colSums(!is.na(rob[-1])))

Analysable resamples per arm:
M1 M2 M3 M4 M5 M6 M7 M8 
30 30 30 30 30 30 30 30 

line("3.7b FACTORIAL MODEL")

========================================================================
3.7b FACTORIAL MODEL
========================================================================

print(anova(fm, type=3))

Type III Analysis of Variance Table with Satterthwaite's method
                                Sum Sq  Mean Sq NumDF DenDF  F value    Pr(>F)    
topology                      0.018368 0.018368     1   203  29.1913 1.827e-07 ***
imputation                    0.000068 0.000068     1   203   0.1078    0.7430    
balancing                     0.112374 0.112374     1   203 178.5938 < 2.2e-16 ***
topology:imputation           0.000000 0.000000     1   203   0.0000    0.9946    
topology:balancing            0.050988 0.050988     1   203  81.0342 < 2.2e-16 ***
imputation:balancing          0.000005 0.000005     1   203   0.0077    0.9303    
topology:imputation:balancing 0.000097 0.000097     1   203   0.1539    0.6952    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

line("3.8  CALIBRATION")

========================================================================
3.8  CALIBRATION
========================================================================

print(calib)

                           Model      Correction  Brier Brier_ref     BSS Calib_intercept Calib_slope
1                   M1: Baseline            none 0.0504    0.0465 -0.0824          -0.206       0.257
2          M2: SICE + Discretize            none 0.0509    0.0465 -0.0941          -0.217       0.240
3           M3: Continuous Nodes            none 0.0561    0.0465 -0.2049          -0.643       0.344
4          M4: SMOTE on Baseline            none 0.1735    0.0465 -2.7279          -3.012       0.153
5          M4: SMOTE on Baseline prior-corrected 0.0489    0.0465 -0.0519          -0.041       0.174
6  M5: ALL (SICE + Cont + SMOTE)            none 0.1756    0.0465 -2.7739          -3.282       0.174
7  M5: ALL (SICE + Cont + SMOTE) prior-corrected 0.0542    0.0465 -0.1651          -0.311       0.193
8          M6: SICE + Continuous            none 0.0564    0.0465 -0.2132          -0.651       0.345
9         M7: Continuous + SMOTE            none 0.1768    0.0465 -2.8007          -3.401       0.165
10        M7: Continuous + SMOTE prior-corrected 0.0552    0.0465 -0.1859          -0.430       0.188
11              M8: SICE + SMOTE            none 0.1750    0.0465 -2.7611          -3.049       0.139
12              M8: SICE + SMOTE prior-corrected 0.0493    0.0465 -0.0593          -0.077       0.169

cat(sprintf("\nPrevalence-only reference Brier = %.4f\n", mean((PI_POP-as.numeric(as.character(testset$stroke)))^2)))

Prevalence-only reference Brier = 0.0465

line("3.9  BENCHMARK")

========================================================================
3.9  BENCHMARK
========================================================================

print(bench)

                                    Model Test_AUC
1 Bayesian network (M3: Continuous Nodes)   0.8109
2       Logistic regression (transparent)   0.8188
3                        XGBoost (opaque)   0.8222

cat(sprintf("XGBoost best_iteration = %d; scale_pos_weight = %.3f\n", xgb$best_iteration, SPW))
cat("\nTop SHAP features:\n"); print(round(head(imp,8),5))

Top SHAP features:
                 age                  bmi    avg_glucose_level         genderFemale        hypertension1 
             0.85756              0.15479              0.10572              0.02464              0.02349 
   work_typeGovt_job smoking_statussmokes          (Intercept) 
             0.01622              0.01375              0.01166 

line("END OF RESULTS PACK")

========================================================================
END OF RESULTS PACK
========================================================================