The dataset comprised 7 variables for prediction, obtained from clinical interviews, physical examination, and routine enzymatic blood laboratory assays. These data were recorded on standard clinical practice by physicians specialized in MODY management during patient assessment and follow-up. These included two binary features (sex and family history of diabetes) and five numerical variables (body mass index —BMI—, blood glycosylated hemoglobin —HbA1c— levels, age at diagnosis, fasting blood glucose levels, and post-load blood glucose levels 120 min after the consumption of 75 g of glucose). These variables are routinely used in the medical practice of Hospital de Clínicas José de San Martín to distinguish persons with MODY from persons with other diabetes phenotypes or from persons without diabetes, and to partially guide medical suspicion regarding which persons are candidates for specific pathogenic variants in MODY-related genes. The selection of these 7 variables was based on prior clinical knowledge of their association with certain cases of GCK-MODY and HNF1A-MODY, as well as on their availability in the hospital records.

Given that MODY variants are not mutually exclusive, binary classification models were trained separately for each target class. Therefore, the dataset was partitioned based on Sanger sequencing results to create distinct subsets for GCK-MODY (162 positive and 60 negative) and HNF1A-MODY (48 positive and 74 negative).

The assessed algorithms can be broadly categorized into traditional parametric models and more flexible non-parametric approaches. The parametric models assume an explicit functional form linking the input features to the class membership probability. This group includes Logistic Regression, which models the probability of the positive class through a sigmoid transformation of a linear predictor; the Perceptron and the Stochastic Gradient Descent Classifier (SGD), which iteratively optimize linear decision functions by minimizing a specified loss function; Linear Discriminant Analysis (LDA), which derives linear decision boundaries under the assumption of Gaussian class-conditional distributions with shared covariance matrices; and Partial Least Squares Discriminant Analysis (PLS-DA), which projects the predictors onto a latent space that maximizes covariance with the response variable. In contrast, non-parametric or semi-parametric models impose fewer assumptions on the data-generating process and are capable of capturing complex non-linear relationships. These include Support Vector Machine (SVC) with kernel functions, which construct optimal separating hyperplanes in high-dimensional feature spaces; K-Nearest Neighbors (KNeighbors), which performs classification based on local similarity among observations; Gaussian Naive Bayes (GaussianNB), which combines simple probabilistic assumptions with empirical estimates of feature distributions; and ensemble-based methods such as Random Forest and XGBClassifier (XGBoost), which aggregate multiple decision trees to enhance predictive accuracy and generalization. Together, these 10 complementary modeling paradigms enabled a comprehensive assessment of classification performance across a wide range of statistical assumptions and learning capacities.

For each machine learning algorithm, a hyperparameter search space was defined (Supplementary Table S1) and systematically explored using grid search implemented via Ray Tune in combination with GridSearchCV. The model configuration achieving the highest Receiver Operating Characteristic Area Under the Curve (ROC AUC) on the training data was identified as the best-performing model. This model was then evaluated on the validation split to compute performance metrics, including ROC AUC, sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and Brier score, which were retained for subsequent analysis. Model validation was performed using bootstrap resampling to generate multiple randomized training and validation splits, enhancing robustness against data variability. For each bootstrap iteration, a hyperparameter search was conducted. Finally, the metrics were averaged across bootstrap iterations to obtain stable and unbiased performance estimates.

Univariate analyses of the zero-imputation GCK-MODY dataset are summarized in Table 1A and provide an initial characterization of differences between diagnostic groups. Three metabolic variables—fasting glucose, post-load glucose, and HbA1c—were significantly associated with diagnostic status, with individuals carrying a positive GCK-MODY diagnosis exhibiting higher mean values across all three parameters. These findings establish a baseline description of the variables most strongly linked to diagnostic outcomes and motivate further exploratory analyses aimed at understanding how these differences manifest across the data.

To visually assess these univariate associations, boxplots (Supplementary Figures S1A–S2A) were first examined, revealing observable shifts in central tendency between diagnostic classes for several continuous predictors in the Full dataset. These patterns were complemented by distribution plots (Supplementary Figure S3A), which further illustrated differences in the overall shape and spread of both continuous and binary variables across groups. Despite these visible shifts, substantial overlap between distributions was consistently observed, indicating that individual variables, when considered in isolation, provide limited discriminative power for separating diagnostic classes.

Building on these observations, correlation maps (Supplementary Figure S4A) were used to examine relationships among predictors and their association with diagnostic status. Age at diagnosis showed negative correlations with most variables, whereas family history of diabetes, HbA1c, fasting glucose, and post-load glucose were positively correlated with a positive GCK-MODY diagnosis, reinforcing the trends identified in the univariate analyses. Finally, principal component analysis biplots (Figure 1A) provided a multivariate overview of the data structure. While some tendency toward clustering by diagnostic status was apparent, considerable intra-class dispersion and partial overlap between groups persisted. Collectively, these results highlight a complex and partially overlapping feature space, underscoring the need for multivariate and machine learning approaches to effectively address the classification task.

A total of 2,500 machine learning models were trained and evaluated for GCK-MODY prediction. Performance metrics are summarized in Table 2A for ROC AUC and in Supplementary Tables S2A–S6A for sensitivity, specificity, PPV, NPV, and Brier score. The final hyperparameter configurations for the best-performing models are reported in Supplementary Table S7.

GaussianNB classifier achieved the highest predictive performance, with a ROC AUC of 0.724, using the dataset processed with full MICE imputation and oversampling. This result reflects moderate discrimination ability. KNeighbors, LDA, logistic regression, and PLS-DA models achieved comparable performance. In contrast, the XGBoost yielded the lowest predictive power, with a ROC AUC of 0.658 under the same preprocessing conditions.

Across models, oversampled datasets consistently outperformed their unbalanced counterparts, and MICE-based imputations generally resulted in better performance than KNN imputations. Calibration analysis based on the Brier score showed that GaussianNB exhibited improved calibration when trained on oversampled data, whereas oversampling led to increased miscalibration for all other models.

The best-performing GCK-MODY model achieved a sensitivity of 0.885 and a specificity of 0.570, with a PPV of 0.853 and an NPV of 0.656. These results indicate good ability to correctly identify true GCK-MODY cases, albeit at the expense of limited performance in correctly identifying true negatives.

Overall, these results indicate that machine learning models trained on routinely collected clinical variables can discriminate GCK-MODY cases from non-cases better than chance, albeit with moderate accuracy. The observed performance supports the feasibility of predicting GCK-MODY subtypes.

To interpret the predictions of the best-performing GCK-MODY model, SHAP-based explainability analyses were conducted. Figure 2A summarizes the overall contribution of each variable to the model’s predictions. age at diagnosis, fasting glucose and HbA1c emerged as the most influential factors driving positive GCK-MODY predictions. Conversely, fasting glucose, BMI, post-load glucose, and age at diagnosis were the primary contributors to negative predictions.

Although the relative importance of predictors varied across individual samples, the model consistently relied on variables that were also statistically significant in univariate analyses, namely fasting glucose, post-load glucose, and HbA1c. This concordance between traditional statistical analysis and model explainability supports the internal coherence of the predictive framework. Additionally, Supplementary Figure S5 illustrates SHAP explanations at the individual level, demonstrating how specific variables influenced model predictions for selected samples and enabling transparent case-by-case assessment.

Importantly, the alignment between SHAP-derived feature contributions and clinically established biomarkers of GCK-MODY indicates that the model’s predictions are driven by meaningful and interpretable clinical signals. This supports the hypothesis that explainable machine learning approaches can provide transparent insights into the variables underlying subtype-specific MODY predictions.

Univariate analysis results for the zero-imputation HNF1A-MODY dataset are reported in Table 1B and reveal a markedly different pattern compared with the GCK-MODY cohort. None of the evaluated variables reached statistical significance at the univariate level. Nonetheless, post-load glucose and age at diagnosis exhibited trends toward association with diagnostic status, with p-values of 0.063 and 0.132, respectively, suggesting weak and inconclusive signals that warrant further multivariate investigation.

Univariate visualizations were subsequently examined to further characterize these patterns. Box plots of continuous variables and summaries of binary predictors (Supplementary Figures S1B–S2B) showed only modest differences in central tendency between diagnostic classes, with substantial overlap across groups. Distribution plots (Supplementary Figure S3B) reinforced this observation, indicating limited discriminative capacity when individual variables were considered in isolation. Consistent with these findings, correlation maps (Supplementary Figure S4B) revealed a structure broadly similar to that observed in the GCK-MODY analysis, with age at diagnosis negatively correlated with most predictors. In contrast to the GCK-MODY dataset, however, post-load glucose emerged as the only variable exhibiting a noticeable positive correlation with HNF1A-MODY diagnosis, underscoring a narrower set of potentially informative features.

Principal component analysis biplots (Figure 1B) further corroborated the limited separability suggested by the univariate analyses, showing pronounced overlap between diagnostic classes and an absence of well-defined clustering by disease status. Collectively, these descriptive, visual, and exploratory findings highlight the greater difficulty of class separation in the HNF1A-MODY dataset and reinforce the need for multivariate modeling approaches to capture subtle and potentially nonlinear relationships within the data.

For HNF1A-MODY, 2,500 models were similarly trained and evaluated. Performance metrics are summarized in Table 2B and Supplementary Tables S2B–S6B. The Random Forest classifier achieved the highest ROC AUC of 0.712, corresponding to moderate classification performance. This result was obtained using the full MICE imputation with unbalanced classes, although comparable performance was observed across datasets with full imputations.

KNeighbors models achieved their best results on datasets with half imputations, while the Perceptron classifier showed the lowest overall performance, reaching a maximum ROC AUC of 0.555 on the zero-imputation oversampled dataset. These findings underscore the sensitivity of model performance to imputation strategy and class balancing. Calibration analysis revealed evidence of miscalibration across all models, which generally worsened with oversampling.

The best-performing HNF1A-MODY model achieved a sensitivity of 0.468 and a specificity of 0.757, with PPV and NPV values of 0.579 and 0.682, respectively. These metrics indicate limited ability to identify true positive cases, but a moderate capacity to correctly rule out HNF1A-MODY in negative individuals.

Overall, these findings indicate that machine learning models trained on routinely collected clinical data can discriminate HNF1A-MODY cases from non-cases better than chance, although with lower sensitivity compared to GCK-MODY. This supports the feasibility of the proposed modular framework while highlighting subtype-specific differences in achievable predictive performance.

SHAP-based explainability analysis was applied to the best-performing HNF1A-MODY model to elucidate the variables driving its predictions. As shown in Figure 2B, age at diagnosis, fasting glucose, post-load glucose, and HbA1c were the most influential variables contributing to negative predictions. In contrast, BMI, post-load glucose, and age at diagnosis had the greatest impact on positive predictions.

Even though univariate analyses did not identify statistically significant predictors for HNF1A-MODY, SHAP results indicate that the model consistently assigned substantial importance to post-load glucose and age at diagnosis. Notably, these variables exhibited trends toward statistical significance in descriptive analyses, suggesting that multivariate modeling was able to capture subtle patterns not evident in univariate testing alone.

Despite the absence of statistically significant predictors in univariate analyses, SHAP-based explanations suggest that the model relied on clinically plausible variables to inform its predictions. This indicates that explainable machine learning can uncover multivariate patterns relevant to HNF1A-MODY while maintaining transparency in how individual clinical features contribute to subtype-specific predictions.