Bayesian Hierarchical Modeling for Predicting Treatment Outcomes in Breast Cancer Patients

Breast cancer is the most prevalent malignancy worldwide with approximately 2.3 million new cases annually. In Kenya, it disproportionately affects younger women (35-50 years) and represents the leading cancer diagnosis. Current prediction models inadequately quantify uncertainty in tumor responses, leading to suboptimal clinical decision-making due to insufficient integration of clinical variables and failure to account for biological variability and institutional heterogeneity. This study developed a Bayesian hierarchical modeling framework to predict pathological complete response (pCR) to neoadjuvant chemotherapy in breast cancer patients. Specific objectives were to develop a Bayesian hierarchical framework incorporating clinical, pathological and treatment variables, compare performance with existing frequentist models, evaluate robustness through sensitivity analyses and obtain probability distributions of treatment outcomes. We conducted retrospective analysis of de-identified electronic health records from 284 patients across 12 Kenyan treatment centers (2018-2022). A Bayesian two-level hierarchical logistic regression model incorporated tumor stage, hormone receptor status, HER₂ expression and treatment center random effects. Markov Chain Monte Carlo methods with Gibbs sampling estimated posterior distributions using OpenBUGS software. Three independent chains with 20,000 iterations each, 12,000 burn-in periods, and thinning intervals of 15 ensured convergence. Multiple imputation addressed missing data patterns. The study population had median age 52 years with 83.4% presenting Stage II – III disease and 38.0% achieving pCR. The enhanced Bayesian hierarchical model demonstrated superior predictive accuracy (AUC = 0.837, Brier score = 0.167) compared to frequentist approaches. Significant prognostic effects were identified for tumor stage (Stage IV vs I: OR = 24.38, 95% CI: 6.60-93.11), hormone receptor positivity (OR = 0.31, 95% CI: 0.15-0.66), and HER₂ positivity (OR = 2.33, 95% CI: 1.04 – 5.32). Treatment center heterogeneity accounted for 9.4% of outcome variability after covariate adjustment. Comprehensive sensitivity analysis demonstrated 97.0% overall stability across validation scenarios. Molecular subtype analysis revealed triple-negative tumors showed 62.4% pCR probability while HR+/HER₂– tumors showed 18.7% probability. The model enables risk-stratified treatment planning with clinically meaningful probability thresholds, optimized resource allocation through center-specific monitoring, and personalized decision-making in resource-constrained settings. The exceptional robustness supports clinical implementation for improving treatment outcomes through enhanced uncertainty quantification and evidence-based treatment selection.