Exploratory data analysis
EDA was performed using matplotlib, pandas, and seaborn libraries to understand the underlying distributions and correlations in the dataset. The dataset has 917,049 rows with information on resistance interpretations of MIC values when tested against a panel of 50 antibiotic drugs for 324 bacterial pathogen species collected across 83 countries from 2004 to 2022, along with patient demographics. A subset of the dataset with 589,998 isolates contains data on the presence or absence of genotype markers such as CTXM, predominantly for ESBL isolates of the Enterobacteriaceae family.
The dataset contained the highest frequency of data collected from the United States of America with the majority share (31%), followed by Spain (12%), France (11.5%), Germany (10.9%), and China (6.1%), aligning with the notion that better monitoring resources are available in economically developed countries17. AMR is a looming threat in Sub-Saharan Africa18; however, there is a significant underrepresentation of data from this region. This disparity highlights the need for increased investments in AMR surveillance infrastructure in low- and middle-income countries (LMICs). Shown in Supplementary Figure S1 is the overall global distribution of antibiotic-resistant bacteria captured in the Pfizer ATLAS Antibiotics dataset.
Missing data is a major problem during prediction; a heatmap (shown in Fig. 1) was generated to understand the number of samples in each column where no data was available. In the heatmap, yellow indicates the absence of data, while purple indicates the presence. No data was missing for the initial columns, such as Isolate ID, Family, and Gender. In contrast, a major part of the data was absent for genetic markers like AMPC, CTXM1, etc., which is reasonable given the extensive resources required for genetic marker analysis in routine surveillance systems. This visualization highlights the need for methods to handle missing data while building machine learning models for practical applications. Imputation is a powerful technique for effectively handling missing data in machine learning, as used by Mutisya et al. 202415 in their study to predict AMR from surveillance data. Although it increases the prediction accuracy of the models, from a clinical perspective, it needs to be done with caution as it could mislead decision-making processes if assessment of the imputation methods is not done19.
Missing data patterns in the dataset for each column are shown on the horizontal axis. Yellow represents missing data, and black represents available data.
Temporal analysis was performed to further understand the year-wise distribution of total samples and resistant samples between the years 2004 and 2022. The total number of samples showed a significant increase up to 2014, reflecting better surveillance and data collection efforts over time. However, a minor decline in total samples can be observed after 2016. The resistance rate can be deduced visually by comparing the two-line plots (Fig. 2). If the resistant line grows at a similar rate to the total line, it implies that the proportion of resistant cases remains largely unchanged, indicating persistent resistance. If the gap between the total and resistant samples increases, it may indicate effective interventions or reduced resistance levels. In contrast, if the lines converge, it could signal an increase in resistance rates. The trend in resistant samples mainly follows the trend of total samples, suggesting a persistent antimicrobial resistance. This could also occur due to sampling bias.
Year-wise Total Samples vs. Resistant Samples from 2004 to 2022.
Correlation analysis was conducted on the full dataset, including genotypic features indicating the presence or absence of β-lactamase genes such as TEM, AMPC, and NDM20,21. These genes are particularly relevant in ESBL-producing organisms like Escherichia coli and Klebsiella spp., where CTX-M, SHV, and TEM serve as hallmark markers22. Co-occurrence of AmpC and carbapenemase genes (e.g., NDM, OXA, IMP) suggests extensive multidrug resistance, limiting therapeutic options23.
To examine associations between bacterial, demographic, genotypic, and phenotypic variables, a Cramér’s V correlation heatmap was generated (Fig. 3). This method, suited for categorical data, revealed strong correlations between species and family, and moderate-to-strong associations between key resistance genes and bacterial taxa (e.g., CTX-M1, SHV, TEM). Carbapenemase genes such as NDM and VIM also showed moderate associations with specific species, indicating species-driven gene dissemination patterns24.
Antibiotic susceptibility patterns were moderately correlated with AMR outcomes, validating the phenotypic influence of tested antibiotics. In contrast, demographic features like age group, gender, and country exhibited weak correlations with genotypic resistance or AMR phenotype. However, the phenotypic resistance class (e.g., MRSA) is predictive of whether the isolate is Susceptible, Intermediate, or Resistant to an antibiotic, with high correlation values. This reflects biological causality, as phenotypes directly influence susceptibility. Source, Species, and Family also show moderate correlations with Phenotype and AMR. This implies that ecological or biological niches (e.g., certain sources like urine, blood, or respiratory) and taxonomic identity affect resistance patterns.
Exploratory analysis of year and country data showed minimal to moderate associations with pathogen and antibiotic variables. This suggests some geographic and temporal variation in testing or prevalence, but not strong enough to drive stratified modeling in this study. These findings support the inclusion of genetic features for improved model interpretability and reinforce the need for global, standardized surveillance of resistance trends.
Cramer’s V correlation plot for all features in the dataset.
Model performance
The dataset was pre-processed and was followed by training and testing classification models for the prediction of antibiotic resistance into Susceptible, Intermediate, and Resistant categories. Phenotype-Only and Phenotype + Genotype datasets were derived from the original dataset and processed individually to avoid selection bias since only a subset of the dataset contains genotype data. The models trained included simplistic models like Logistic regression as well as complex ensemble models like XGBoost. The model performance metrics for the initial set of models with all features included and implemented are given in Table 1 for the Phenotype-only and Phenotype + Genotype datasets. Simpler models like logistic regression and KNN have lower accuracies, whereas more complex models like XGBoost and Gradient Boost have better prediction ability and higher scores on all metrics. For the Phenotype-only dataset, SVM has the highest precision. For both datasets, XGBoost was the best-performing model, with its AUC being the highest (95%) and its accuracy slightly enhanced with hyperparameter tuning and ten-fold cross-validation. Supplementary Table S2 shows the parameters used for model development and the tuning parameters used for hyperparameter tuning.
Shown in Fig. 4 are the ROC-AUC curves of the XGBoost models. ROC-AUC measures the model’s ability to distinguish between classes. A higher AUC indicates a better-performing model, with 1.0 being the best possible score and 0.5 representing random guessing. To handle the multi-class nature of the problem, the One-vs-Rest (OvR) strategy was implemented. For each class, the model was trained to distinguish it from the others, and binary classification was performed on each one. The ROC curve was plotted by calculating the False Positive Rate (FPR) and True Positive Rate (TPR) for each class. The AUC for each class was then calculated, with a higher AUC indicating a better performance at distinguishing that class from the others. To statistically evaluate differences in model performance, the DeLong test was applied to compare the AUC values between models. The results indicated no statistically significant difference in AUC values, suggesting that the models exhibit similar discriminatory power. Additionally, the bootstrapped confidence intervals around the AUC estimates were narrow, indicating stability and low variability in the AUC across resampled datasets. While this supports the robustness of the models’ performance, it also highlights the limitation of relying solely on AUC for model selection. Given the lack of significant difference in AUC, additional metrics—such as precision, recall, and F1-score—as well as model interpretability, were considered to enable a more comprehensive comparison of model behavior and utility.
ROC-AUC curves for XGBoost models for (a) Phenotype-Only and (b) Phenotype + Genotype datasets. Classes 0,1, and 2 are Intermediate, Resistant, and Susceptible, respectively.
Data balancing
The accuracy of prediction is lowest for the Intermediate category for both Phenotype-only and Phenotype + Genotype datasets, as seen in the confusion matrices of the XGBoost models (Fig. 5). This is because of data imbalance, with the number of intermediate samples being considerably lower compared to resistant and susceptible samples, which could affect the model’s ability to generalize effectively (S/I/R data distribution shown in Supplementary Figure S2). Imbalanced datasets could cause ML models to bias toward the majority class, misclassify rare but clinically important classes, and underperform in metrics like recall and F1-score for minority classes. Data balancing using under-sampling of susceptible samples and oversampling of intermediate samples done on the XGBoost models improved recall.
Confusion matrices of the XGBoost models for (a) Phenotype-Only and (b) Phenotype + Genotype datasets.
Feature importance
For models like Logistic Regression, Random Forest, Gradient Boosting, AdaBoost, and XGBoost that support feature analysis, the corresponding scores that indicate how much each feature contributes to the model’s performance were extracted. These feature importance scores reflect how often and how effectively a feature contributes to the decision-making process in the model’s trees. These scores were used to determine which groups of features had the most significant impact on the model’s predictions (Fig. 6).
Feature importance plots of XGBoost models for (a) Phenotype-Only and (b) Phenotype + Genotype datasets.
SHAP (SHapley additive exPlanations) analysis
SHAP analysis25 for the XGBoost model trained on the Phenotype-only dataset (summary plot shown in Fig. 7) revealed that individual features such as Phenotype_ESBL, Phenotype_MRSA, and antibiotics like tigecycline, vancomycin, and oxacillin were among the most influential in predicting AMR phenotypes. These features are known resistance indicators and reflect the model’s reliance on clinically relevant patterns26,27,28,29,30. Most top features contributed to the prediction of the intermediate class, while fewer features drove predictions for susceptible, suggesting class-wise differences in feature influence. The SHAP plot enhances interpretability by identifying which phenotypic or treatment features the model used most heavily, supporting transparency and clinical relevance.
The SHAP summary plot for the Phenotype + Genotype dataset (Fig. 7) provides a global view of the most influential features driving model predictions for antimicrobial resistance phenotypes—Susceptible (Class 2), Resistant (Class 1), and Intermediate (Class 0). Features such as Antibiotic_Colistin, Antibiotic_Ampicillin, and Antibiotic_Tigecycline exhibited the highest average SHAP values, indicating strong contributions to the model’s predictive accuracy. Notably, Colistin and Ampicillin showed the greatest impact on predicting the Intermediate phenotype (Class 0), as seen by the dominance of green bars31,32. Conversely, features like Ceftazidime avibactam, Cefaroline, and resistance genes such as OXA_NEG and NDM_NEG contributed significantly to identifying Resistant phenotypes (Class 1), suggesting the model correctly leveraged key molecular determinants of resistance24,33. Susceptible predictions (Class 2) were more associated with features such as Escherichia coli, Cefepime, and absence of resistance genes (KPC_NEG, TEM_NEG), indicating that the model identified typical patterns associated with non-resistant strains34,35. These findings not only validate the model’s behavior but also enhance interpretability by aligning with known biological mechanisms. Such model explanations are essential in AMR contexts where understanding the why behind a prediction is as important as the prediction itself. Although genotype features offered modest overall performance gains, they showed disproportionately high importance in ESBL and carbapenem resistance cases. This underscores the clinical value of including molecular features in high-priority resistance scenarios.
SHAP plots of XGBoost models for (a) Phenotype-Only and (b) Phenotype + Genotype datasets, where classes 0, 1 and 2 are Intermediate, Resistant and Susceptible, respectively. SHAP summary plot illustrating the contribution of each feature to the model’s predictions. Features are ranked by importance (y-axis), and the horizontal distribution shows their impact on the prediction output.
