This study systematically investigated and compared eight ML models for predicting VIN in pediatric populations. Three principal findings emerge from our analysis. First, XGBoost, AdaBoost, RF, and LightGBM achieved comparable high predictive performance with AUROCs of 0.88–0.90, with no statistically significant differences among these four models (all P > 0.05 by DeLong's test). XGBoost was selected as the final model based on its favorable balance of discrimination (AUROC = 0.90), precision-recall performance (AUPR = 0.74), and calibration (Brier score = 0.07), outperforming traditional logistic regression approaches reported in prior pediatric studies (AUROC 0.66–0.83) [4,5,6,7]. Second, SHAP-based interpretability analysis identified eGFR, amphotericin B co-administration, and ETC as the most influential predictors, with eGFR alone contributing a mean importance value of 0.204. It is important to note that while baseline eGFR did not differ significantly between VIN and non-VIN groups in univariate analysis (P = 0.15), XGBoost exploited the full continuous distribution of eGFR and its non-linear interactions with other variables, enabling superior risk stratification beyond simple group comparisons. Third, a parsimonious six-feature model maintained comparable predictive accuracy (AUROC 0.87, P > 0.05 vs. full model by DeLong's test) with acceptable calibration (Brier score 0.08, ECE 3.8%). While decision curve analysis suggested theoretical net benefit under modeling assumptions, this does not establish clinical effectiveness or readiness for deployment. External validation, prospective implementation studies, and interventional trials are required to evaluate whether model-guided risk stratification improves patient outcomes, resource utilization, or nephrotoxicity rates in real-world practice.

Strengths. Our study possesses several methodological strengths. The large sample size (1,452 children) represents one of the largest pediatric VIN machine learning cohorts reported to date, providing robust statistical power for model development and internal validation. The systematic comparison of eight diverse machine learning algorithms, including tree-based methods (XGBoost, LightGBM, RF, AdaBoost, DT), instance-based methods (KNN), NNET, and SVM, ensures that the selected model represents a data-driven choice rather than an arbitrary choice. The application of SMOTE to address class imbalance (14.2% VIN rate) reduces bias toward the majority class, a critical consideration often neglected in adverse event prediction studies. The use of 100 bootstrap replicates for feature importance estimation enhances the stability and reliability of predictor rankings. Most importantly, the integration of SHAP analysis at both the population summary level (beeswarm plot) and individual patient level (waterfall plots) partially overcomes the "black-box" limitation that has historically impeded clinical adoption of machine learning models [15]. The derivation of a simplified six-feature model via stepwise reduction guided by importance ranking, with non-inferiority confirmed by DeLong's test, balances predictive accuracy with clinical practicality, a framework for future prospective validation rather than established clinical utility. Additionally, we employed rigorous dual-validation methodologies beyond the standard single train-test split. Repeated nested cross-validation with stratified sampling and random-search hyperparameter tuning provided 15 independent performance estimates with full uncertainty quantification, while bootstrap optimism correction offered a conservative, bias-adjusted performance estimate. The remarkable consistency between the primary test-set AUROC (0.872) and nested CV mean (0.869) indicates that the single split did not introduce meaningful optimism bias, likely due to the large sample size (n = 1,452) and stratified randomization. The bootstrap-corrected estimate (0.816) provides clinicians with a realistic lower-bound for expected performance in practice.

Weaknesses. Several limitations of this study should be acknowledged. First, the retrospective single-center design introduces generalizability limitations. Center-specific prescribing practices, including vancomycin dosing protocols, concomitant medication preferences, and nephrotoxicity monitoring frequencies, may limit external validity. The model was developed in a single tertiary children's hospital in China, where weight-based dosing (median IVDW 10.2 mg/kg/dose) is standard practice; institutions with alternative strategies (e.g., area-under-the-curve-guided dosing) may experience different nephrotoxicity risk profiles. The high prevalence of amphotericin B co-administration in our cohort (66/1,452, 4.5%) reflects local antifungal prescribing patterns, and hospitals with different formulary preferences or antifungal stewardship programs may exhibit different predictor-outcome relationships. Of the initial cohort, approximately one-third of patients were excluded due to missing height data, which precluded eGFR calculation. While the missing height data mechanism appears primarily administrative rather than clinically driven, the substantial exclusion rate may introduce potential selection bias and limit the generalizability of our findings to populations with complete anthropometric documentation. External validation in independent pediatric cohorts across diverse geographic, ethnic, and clinical settings, ideally involving prospective multi-center validation and integration of real-time monitoring data, is therefore essential prior to any clinical application. Second, we did not adjust for illness severity scores (e.g., PRISM) due to incomplete data capture in our EHR system. Consequently, disease severity may confound the association between predictors (e.g., eGFR, concomitant medications) and VIN risk, as sicker patients are more likely to receive nephrotoxic agents and develop renal dysfunction. This unmeasured confounding represents a significant limitation: the observed associations between immunosuppressant use, amphotericin B exposure, and VIN risk may partially reflect underlying disease severity rather than direct causal effects of these medications. Future studies should incorporate comprehensive severity scoring systems to disentangle medication effects from illness severity. Third, the final model demonstrated only moderate sensitivity (0.64), indicating that a substantial proportion of VIN cases may go undetected, thus, ongoing clinical vigilance remains necessary even when the model predicts a low risk. Fourth, while SHAP analysis provides post-hoc explanations of model behavior, it cannot establish causal relationships; Accordingly, the identified associations should be interpreted as predictive rather than etiological. Fifth, our definition of nephrotoxicity based on serum creatinine changes may underestimate true kidney injury, particularly in young children where muscle mass and creatinine generation rates vary substantially with age. The Schwartz formula used for eGFR estimation partially accounts for age through height, but does not fully capture developmental changes in renal function. Consequently, the fixed threshold of ≥ 0.5 mg/dL or ≥ 50% increase may be less sensitive in younger children with lower baseline creatinine, potentially leading to misclassification of mild kidney injury. The absence of newer renal injury biomarkers, such as cystatin C, further reduces diagnostic certainty and presents an opportunity for future model refinement [27]. Sixth, features with missing data exceeding 25% were excluded during model development, which may have led to the overestimation of importance of retained predictors and introduced potential selection bias. Seventh, we did not perform temporal validation, which evaluates model performance on data collected after the training period and is particularly valuable for detecting dataset shift or temporal drift in clinical practice patterns. Our 7-year study period (December 2017–November 2024) encompassed evolving vancomycin dosing protocols, changing nephrotoxicity surveillance practices, and the COVID-19 pandemic period (2020–2022), all of which could theoretically introduce temporal heterogeneity. Eighth, the calibration slope of 0.661 indicates moderate under-confidence, where predicted probabilities are shrunk toward the baseline event rate (14.2%). While this conservatism does not impair the model's primary objective of identifying high-risk patients (predicted probabilities ≥ 0.142 remain highly specific for true VIN), it may lead to under-estimation of absolute risk for individual patients. A calibration slope below 1 implies systematic underestimation of predicted probabilities in high-risk patients, which is of particular concern for a nephrotoxicity risk stratification tool where under-triage of the highest-risk cases carries significant clinical consequences. Additionally, the lowest-risk bin showed the largest calibration error (MCE 8.5%), attributable to the imbalanced class distribution; clinicians should interpret very low predicted probabilities (< 0.05) as indicating "baseline risk" rather than absence of risk. Ninth, the model requires ETC measured within 48 h, limiting its utility for very early risk stratification at the time of vancomycin initiation. A pre-trough model excluding ETC would require separate development and validation.

Among the top-performing models, XGBoost, LightGBM, AdaBoost, and RF achieved comparable AUROCs (0.88–0.90) with no statistically significant differences (all P > 0.05 by DeLong's test). XGBoost and LightGBM demonstrated identical discrimination (AUROC = 0.90, AUPR = 0.74) and near-identical performance across all metrics (differences < 0.03), indicating practical equivalence in predictive capacity.

The superior performance of tree-based ensembles (XGBoost, RF, LightGBM, AdaBoost; AUROCs 0.88–0.90) compared to SVM (0.75), KNN (0.80), and NNET (0.80) likely reflects the presence of non-linear relationships and feature interactions in this dataset, which tree models capture through recursive partitioning. The extremely low sensitivity of SVM (0.02) suggests that the radial kernel and hyperparameter configuration were poorly suited for this imbalanced dataset, resulting in a model that almost never predicted VIN. NNET's underperformance (F1 0.46) may reflect insufficient sample size for robust neural network training or suboptimal architecture selection given our threefold cross-validation constraint.

In this imbalanced cohort (14.2% VIN), AUROC can be misleadingly optimistic due to the large true-negative population. AUPR provides more informative assessment of positive-class prediction quality. XGBoost and LightGBM achieved equivalent AUPR (0.74), but LightGBM showed higher sensitivity (0.61 vs 0.58) and accuracy (0.92 vs 0.91). We selected XGBoost based on its superior calibration (Brier 0.07 vs 0.08) and well-established interpretability tools, though this choice involved trade-offs in sensitivity. Additionally, XGBoost's mature regularization framework and extensive validation in adverse drug reaction prediction literature support its reliability for clinical implementation [13, 28,29,30].

Feature importance and SHAP analyses further identified eGFR, amphotericin B co-administration, and ETC as the features with largest predictive contribution. It is critical to emphasize that SHAP values explain model behavior, not causal mechanisms. For example, amphotericin B co-administration may be a proxy for severe fungal infection rather than a direct nephrotoxic cause; eGFR may reflect underlying renal reserve rather than a modifiable target. Any intervention targeting these predictors requires prospective validation in randomized trials, as SHAP identifies predictive rather than causal relationships. The dominant predictive role of eGFR is consistent with established pathophysiological mechanisms: diminished renal clearance is associated with drug accumulation in pathophysiological studies [3, 31, 32]. However, SHAP identifies predictive rather than causal relationships. Notably, while univariate comparison showed no significant difference in median eGFR between groups (P = 0.15), the XGBoost model captured the predictive signal across the continuous eGFR distribution. SHAP quantification further illustrated this association at the individual level: impaired eGFR showed the largest positive SHAP contribution to predicted risk. This demonstrates that ML models can detect clinically meaningful patterns that traditional univariate analyses may miss. Higher ETC was correlated with increased SHAP values, consistent with prior observations of an association between nephrotoxicity risk and high trough concentration. A previous multicenter trial reported that trough concentrations > 15 mg/L were associated with a threefold higher nephrotoxicity risk in adult [33], while, a meta-analysis confirmed that trough levels ≥ 15 mg/L increased nephrotoxicity by 2.7-fold in pediatric patients [34]. Concomitant administration of nephrotoxic agents, particularly amphotericin B, showed a strong predictive association, consistent with prior findings that combined nephrotoxin exposure substantially elevates VIN susceptibility [35].

The higher sensitivity of the 6-feature model (0.64) compared to the full 34-feature model (0.58) is counterintuitive but consistent with the bias-variance trade-off. We hypothesize that the 28 excluded low-importance features introduced noise that suppressed true positive predictions in the full model, particularly at the operating threshold. This phenomenon aligns with the 'curse of dimensionality' in high-dimensional, modest-sample settings and supports the value of parsimonious feature selection in clinical prediction modeling.

Currently, no consensus or guideline defines the optimal features quantity for clinical prediction modeling. Although expanding feature sets may theoretically improve predictive power by incorporating more information, excessive redundant variables introduce noise, weaken clinical feasibility, and may impair model performance [16]. To enhance clinical practicality, we established a parsimonious six-feature XGBoost model via stepwise selection guided by feature importance, including eGFR, amphotericin B co-administration, ETC, ALB, immunosuppressant use, and acyclovir/ganciclovir co-administration. The simplified model yielded an AUC of 0.87, which was comparable to the full 34-feature model (AUC = 0.90; P > 0.05 by DeLong's test), maintaining satisfactory discrimination while greatly improving clinical operability. Our 6-feature XGBoost model demonstrated discrimination comparable to or exceeding prior pediatric VIN prediction tools in this cohort. Traditional logistic regression models with limited predictors reported AUCs ranging from 0.66–0.83 [4,5,6,7], though direct comparison is limited by differences in study design, population, and validation methods. In comparison with the recent study by Yin et al. [8], our work represents a complementary but distinct contribution. Yin et al. focused on predicting drug exposure (trough concentrations) to generate hypotheses for risk stratification and future preventive interventions, a pharmacodynamic and safety objective. Both studies identify renal function indicators (creatinine, BUN, eGFR) as critical predictors, consistent with the central role of kidney function in vancomycin pharmacology. However, our larger sample size, broader age range (1 month to 18 years vs. < 4 years), and comprehensive algorithm comparison extend prior methodological efforts in pediatric VIN prediction. Decision curve analysis, a theoretical framework for evaluating potential clinical utility under specific assumptions, suggested net benefit across risk thresholds (0.05–0.75). This theoretical benefit does not establish clinical effectiveness, safety, or readiness for implementation, but generates hypotheses for future interventional research examining whether model-guided risk stratification could improve therapeutic drug monitoring, dose optimization, and nephrotoxicity prevention in prospective studies.

For the model to be clinically meaningful, high-risk predictions should trigger specific, actionable responses. Based on the identified predictors, we hypothesize that such actions might include: Enhanced renal monitoring: Increased creatinine monitoring frequency (e.g., every 24–48 h vs. every 72 h) for high-risk patients; Pharmacist-led review: Early consultation for vancomycin dose optimization and drug interaction assessment; Nephrotoxin avoidance: Avoidance of additional nephrotoxic agents (particularly amphotericin B) in high-risk patients; Earlier therapeutic drug monitoring: More frequent trough measurements to guide dose adjustment. However, the optimal intervention bundle, its feasibility, and impact on patient outcomes remain unknown and require evaluation in prospective implementation studies.