Interpretable machine learning algorithms for diagnostic prediction of diabetic retinopathy

Abstract

Background

Diabetic Retinopathy (DR) remains a leading cause of blindness among diabetic patients worldwide, necessitating early and accurate diagnostic interventions. While traditional screening methods rely heavily on manual ophthalmologic evaluations, recent advancements in machine learning (ML) and deep learning (DL) have opened new avenues for automated, scalable, and interpretable diagnostic tools. However, challenges persist in developing models that are not only high-performing but also transparent enough to gain clinical trust.

Objective

This study introduces a novel, standardized, and interpretable ML framework designed specifically to enhance diagnostic efficiency and accuracy for DR risk prediction. By prioritizing model interpretability alongside predictive performance, our approach aims to bridge the gap between cutting-edge AI technology and clinical applicability.

Methods

We evaluated eleven ML algorithms, optimizing hyperparameters via grid search and five-fold cross-validation to identify top-performing models. A key innovation lies in our dynamic weighted voting ensemble (Voting_soft), which integrates multiple classifiers based on model confidence, thereby leveraging the strengths of diverse algorithms. Model performance was rigorously assessed using accuracy, sensitivity, and area under the curve (AUC) metrics, with ROC and PR curves comparing performance across varying training dataset proportions. Crucially, we employed SHAP (SHapley Additive exPlanations) for interpretability analysis, providing clinicians with actionable insights into feature contributions.

Results

Through LightGBM-based correlation analysis and AUC curve determination, fourteen clinical features were identified as optimal predictors. Notably, the CatBoost model achieved superior performance on a 20% test set, while the Extreme Random Tree model demonstrated robustness on a 30% test set. Our dynamic weighted voting ensemble (Voting_soft) outperformed individual models in terms of AUC across both datasets. SHAP analysis revealed that age, triglycerides, sex, and HDL-C were key predictors of DR prevalence, offering clinically meaningful explanations for model decisions.

Conclusions

This study presents a groundbreaking ML-based DR risk prediction system that excels in both accuracy and interpretability. The integration of SHAP analysis not only enhances model transparency but also empowers clinicians with a deeper understanding of diagnostic decision-making, ultimately improving the precision and efficiency of DR screening. Our dynamic voting ensemble approach sets a new benchmark for interpretable, multi-model integration in medical diagnostics.

Keywords

machine learning SHAP diabetic retinopathy optimization prediction

1 Introduction

Diabetic Retinopathy (DR) is one of the most common microvascular complications of diabetes and is the leading cause of new cases of blindness in diabetic patients.¹ Globally, more than 132 million people suffer from DR and it is the leading cause of visual impairment and blindness in the working-age population.² In 2020, more than 103 million people with diabetes were affected by diabetic retinopathy globally, and this number is projected to increase to 160 million by 2045.³ Compared to all other major causes of blindness, diabetic retinopathy is the only disease whose age-standardized prevalence did not decline between 1990 and 2020. Therefore, early detection, accurate diagnosis and effective treatment of DR are important for improving patient prognosis and reducing the burden on society.

With the development of medical technology, particularly the emergence of machine learning (ML)⁴ and deep learning (DL)⁵ techniques, new approaches and methodologies have been introduced for the diagnosis and treatment of DR. Wang Qianwen⁶ constructed a DR risk prediction model based on the XGBoost algorithm and further enhanced its predictive performance using a Stacking ensemble model, achieving an AUC value of 0.924. This study not only validated the potential of ensemble learning algorithms in DR risk prediction but also demonstrated effective strategies for improving prediction accuracy through model fusion. OY JL⁷ explored genetic variables associated with DR using Mendelian randomization and developed prediction models incorporating logistic regression and machine learning algorithms, thereby enhancing the accuracy of DR diagnosis. Sunyoung Kim et al.⁸ utilized electronic medical record data from three universities in South Korea to develop a machine learning algorithm based on extreme gradient boosting for predicting DR risk in patients with type 2 diabetes. The results indicated that the XGBoost model exhibited high accuracy and specificity in both the original and validation datasets. Amanda Luong et al.⁹ proposed a novel retinopathy risk score (RRS) model and validated its effectiveness in identifying low-risk DR patients through comparison with machine learning models (MLMs). Zhao YL et al.¹⁰ integrated computer vision and clinical structured data to construct a multimodal model for predicting DR referrals. The findings demonstrated that the XGBoost-based multimodal model performed exceptionally well in both internal validation and external testing sets, significantly improving the accuracy of DR referral prediction. Li BB et al.¹¹ developed a DR risk prediction model using a self-evolving machine learning approach that combined fundus images and clinical data, and calculated DR risk through a visualization system. Pang YH et al.¹² identified potential biomarkers for early DR, diabetic macular edema (DME), and anti-vascular endothelial growth factor (anti-VEGF) treatment responses using multi-omics data (including lipidomics and metabolomics) and machine learning algorithms. Lukashevich Marina M¹³ classified fundus images of diabetic retinopathy using ensemble learning methods and provided references for developing personalized treatment plans by analyzing the performance of different algorithms.

Despite remarkable advancements in the application of machine learning for diabetic retinopathy, the poor interpretability of machine learning models poses a significant challenge. Clinicians often struggle to comprehend the decision-making processes of these models, which has consequently limited their clinical acceptance and adoption.¹⁴

The goal of this paper is to develop a generic, standardized, and effective interpretable machine learning system for predicting the risk of DR. This study also provides the following contributions: the data were preprocessed using preprocessing techniques such as data cleaning and normalization. LightGBM and Random Forest were used as feature selection methods, respectively, and the number of predictors of DR was determined by combining feature selection methods with three models. Highly applicable machine learning algorithms including Random Forest, Decision Tree, Extreme Gradient Boosting, K Nearest Neighbors, Extreme Random Tree, Support Vector Machines, and Multi-Layer Perceptual Machines are used. Fusion machine learning algorithms are explored to achieve voting integration of classification multi-models with VotingClassifier. Hyperparameters were optimized by grid search method and experiments were conducted using a five-fold cross-validation approach. Evaluated model performance using metrics such as accuracy, precision, recall, F1 score, and area under the curve (AUC) using receiver operating characteristic curves (ROC) as well as PR curves. The powerful SHAP technique is applied to interpret and understand the contribution of each predictor in DR prediction. The rest of the paper is organized as follows: section 2 presents the concepts and methods, section 3 presents the results and analysis, and section 4 gives the conclusion and outlook.

2 Concepts and methods

2.1 Machine learning algorithms

Random Forest(RF),¹⁵ proposed by Breiman, is an integrated learning algorithm based on decision trees, which reduces the risk of overfitting by constructing multiple decision trees and integrating their predictions, the core of which lies in the dual randomness of Bootstrap sampling and random selection of features; Decision Tree(DT)¹⁶ serves as a basic classification model, constructing a tree rule through recursive feature partitioning, deriving ID3, C4.5, CART and other classic variants; Extreme Gradient Boosting(XGBoost)¹⁷ is an optimized implementation of gradient boosting decision tree, developed by T.Q. Chen's team, which improves the training efficiency through the second-order Taylor expansion and regularization terms, and excels in structured data modeling; K-nearest-neighbor algorithm(KNN)¹⁸ is based on the idea of local approximation, and selects the K nearest neighbors through the distance metric to make a voting decision, which is suitable for small to medium-sized low-dimensional data; Extreme Random Tree(ET),¹⁹ as a variant of Random Forest, adopts a random threshold splitting strategy, exchanging small deviations for variance reduction; LightGBM,²⁰ proposed by Microsoft, optimizes large-scale data processing through a histogram algorithm and a leaf growth strategy, striking a balance between speed and accuracy; CatBoost(CAT),²¹ developed by Yandex, is designed to design an ordered boosting mechanism for the category features and has a built-in missing value. boosting mechanism, and built-in missing value processing module; AdaBoost(ADA)²² as a classical Boosting algorithm, by dynamically adjusting the sample weights in series of weak classifiers, and ultimately forming a strong predictor; GBM²³ as a gradient boosting framework, through iterative optimization of the loss function to improve the performance of the model, laying the foundation for the subsequent algorithms, such as XGBoost; Support Vector Machines(SVM)²⁴ Based on statistical learning theory, it searches for the maximum interval hyperplane through kernel trick, which has a natural advantage in high-dimensional data classification; Multilayer Perceptron(MLP),²⁵ as a feed-forward neural network, approximates the complex function through multilayer nonlinear mapping, and shows strong expression ability in the field of image, speech, and so on.

2.2 Algorithm evaluation metrics

In this paper, a comprehensive evaluation index system containing more information is adopted. On the basis of the traditional accuracy rate, the recall rate is introduced to measure the model's completeness in recognizing positive cases, i.e., the ratio of the number of correctly determined positive cases to the total number of actual positive cases; the F1 score, as the reconciled mean value of the precision rate and the recall rate, integrally reflects the model's classification effect on the majority class and the minority class, with particular emphasis on the recognition performance of the minority class. Meanwhile, the model's ability to exclude negative cases is assessed by specificity, i.e., the proportion of correctly identified negative cases to the total number of actual negative cases; positive and negative predictive values quantify the reliability of the samples predicted to be positive and negative cases, respectively. In addition, the ROC curve is used to depict the balance of sensitivity and specificity of the classifier under different thresholds, and the larger the AUC value of the AUC, the better the model's ability to distinguish between positive and negative samples; whereas the Precision-Recall curve focuses on the identification performance of positive cases, and the larger the AP value of the AUC, the better the model's classification performance in the data imbalance scenario. The above metrics together build a multi-dimensional evaluation framework, in which the confusion matrix serves as the basis to support the computation and interpretation of each metric through four basic statistics: TP, FN, FP, and TN.²⁶

2.3 Introduction to SHAP

The SHapley Additive exPlanations(SHAP)²⁷ method, rooted in Shapley value theory from game theory, was proposed by economist Lloyd Shapley in 1953 to address fairness in profit distribution within cooperative games. This study conceptualizes the model's prediction process as a cooperative game, where each feature acts as a participant and the model's predictive outcome represents the collective gain. By calculating Shapley values for individual features, SHAP quantifies the marginal contribution of each feature to the model's prediction, thereby elucidating the decision-making process. Features with positive SHAP values contribute to predicting DR occurrence, while those with negative SHAP values support the prediction of non-DR cases.

3 Results and analysis

3.1 Data sources

For the original data of physical examination in a hospital, preprocessing was carried out in this paper, and the main process included (1) Missing values were filled in using the mean filling method for continuous-value data, and plurality was used for dichotomous attribute values, and the descriptive analyses of these attributes are shown in Table 1. (2) Variable mapping was performed: SEX was assigned to 0 for males and 1 for females, and 1 for those with a history of high blood pressure, and 0 for those without.The BMI was set as follows: 0 = less than or equal to 23.9, 1 = 24.0–27.9, and 2 = greater than or equal to 28.0. The final dataset contained 1006 samples totaling 15 attributes.

Table 1.
Descriptive analysis.

mean std min 25% 50% 75% max

age 50.29 16.92 20 37.25 50 60 93

HDLC 1.15 0.31 0.42 0.92 1.12 1.32 2.5

LDLC 2.71 0.85 0.84 2.1 2.68 3.22 7.98

VLDLC 1 0.72 0.14 0.68 0.85 1.09 11.26

TG 1.9 2.42 0.35 0.88 1.34 2.09 45.84

TC 4.86 1.03 2.41 4.2 4.79 5.38 12.61

PR 80.82 12.54 41 72 79 88 135

DBP 76.89 12.76 45 67 76 85 119

BUN 5.56 1.65 2.21 4.45 5.34 6.37 18.64

UA 339.35 84.57 140.8 280.85 333 394 679

CR 64.11 29.34 30 51.25 62 72 799

	mean	std	min	25%	50%	75%	max
age	50.29	16.92	20	37.25	50	60	93
HDLC	1.15	0.31	0.42	0.92	1.12	1.32	2.5
LDLC	2.71	0.85	0.84	2.1	2.68	3.22	7.98
VLDLC	1	0.72	0.14	0.68	0.85	1.09	11.26
TG	1.9	2.42	0.35	0.88	1.34	2.09	45.84
TC	4.86	1.03	2.41	4.2	4.79	5.38	12.61
PR	80.82	12.54	41	72	79	88	135
DBP	76.89	12.76	45	67	76	85	119
BUN	5.56	1.65	2.21	4.45	5.34	6.37	18.64
UA	339.35	84.57	140.8	280.85	333	394	679
CR	64.11	29.34	30	51.25	62	72	799

3.2 Data correlation analysis and feature selection

Figure 1 gives a diagram of the two processes of data preprocessing. Among them, Figure 1A gives a heat map of the correlation coefficient matrix between features. Figure 1B shows the prediction of diabetic retinopathy based on four models selected with different number of features. The AUC in the figure varies with the increase in the number of features and the error bars indicate the 95% confidence intervals. Here Random Forest was used as the feature selection method and screened in combination with three models (Logistic regression, XGBoost and LightGBM). As the number of features increases, the AUC value gradually increases until all the features are included, the increase in AUC tends to level off. Therefore, 14 clinical features such as AGE,TG were selected as inputs to the final model.

Figure 1.

Plot of data correlation analysis and feature selection curves.

3.3 Machine learning algorithm results

3.3.1 Optimal parameter output

Firstly, the dataset is divided into training and test sets, using 70% of the data as the training set and the remaining 30% of the data as the test set, using the algorithm in subsection 2.1 for the binary classification task, setting multiple hyper-parameter options, and searching for the best hyper-parameter combinations through 5-fold cross validation combined with grid search, and finally training the model with the optimal performance, the results are shown in Table 2.

Table 2.
The optimal hyperparameter combinations.

Models Hyperparameter Optimal value

RF n_estimators’: [50, 100, 200], ‘max_depth': [3, 5, 7, None], ‘max_features': [‘sqrt’, ‘log2’, None],'min_samples_split': [2, 5, 10],'min_samples_leaf': [1, 2, 4],'bootstrap': [True, False],'max_samples': [0.8, 1.0] bootstrap': True, ‘max_depth': None, ‘max_features’: ‘sqrt’, max_samples': 0.8, ‘min_samples_leaf': 1, min_samples_split': 5, n_estimators': 200

DT max_depth': [3, 5, 7, None],'min_samples_split': [2, 5, 10],'min_samples_leaf': [1, 2, 4],'max_features': [‘sqrt’, ‘log2’, None],'criterion': [‘gini’, ‘entropy’],'ccp_alpha': [0.0, 0.1, 0.2] ccp_alpha': 0.0, ‘criterion': ‘gini’, ‘max_depth': 5, max_features': ‘sqrt’, min_samples_leaf': 1, min_samples_split': 5

XGBoost n_estimators': [50, 100, 200],'max_depth': [3, 5, 7],'learning_rate': [0.01, 0.1, 0.2],'subsample': [0.8, 1.0],'colsample_bytree': [0.8, 1.0] colsample_bytree': 0.8,'learning_rate': 0.01,'max_depth': 7,'n_estimators': 200,'subsample':0.8

KNN n_neighbors': [3, 5, 7, 9, 11],'weights': [‘uniform’, ‘distance’],'metric': [‘euclidean’, ‘manhattan’, ‘minkowski’],'p': [1, 2] metric': ‘manhattan’, ‘n_neighbors': 9, ‘p': 1, ‘weights': ‘uniform'

ET n_estimators': [50, 100, 200],'max_depth': [3, 5, 7, None],'min_samples_split': [2, 5, 10],'min_samples_leaf': [1, 2, 4],'max_features': [‘sqrt’, ‘log2’, None],'bootstrap': [True, False],'criterion': [‘gini’, ‘entropy’] bootstrap': False, ‘criterion': ‘entropy’, ‘max_depth': None, ‘max_features': None, ‘min_samples_leaf': 1, ‘min_samples_split': 10, ‘n_estimators': 50

LightGBM ‘num_leaves': [31, 50, 100],'learning_rate': [0.01, 0.05, 0.1],'n_estimators': [100, 200, 500],'min_child_samples': [20, 50, 100],'subsample': [0.6, 0.8, 1.0] ‘num_leaves': 31,'learning_rate': 0.01,'n_estimators': 100,'min_child_samples': 20,'subsample': 0.6

CatBoost n_estimators': [50, 100, 200],'max_depth': [3, 5, 7], ‘learning_rate': [0.01, 0.1, 0.2],'subsample': [0.8, 1.0],'colsample_bylevel': [0.8, 1.0],'border_count': [32, 64, 128],'bagging_temperature': [0, 1, 2],'l2_leaf_reg': [0, 1, 3] bagging_temperature': 0, ‘border_count': 128, ‘colsample_bylevel': 1.0, ‘l2_leaf_reg': 0, ‘learning_rate': 0.2, ‘max_depth': 3, ‘n_estimators': 200, ‘subsample': 1.0

AdaBoost n_estimators': [50, 100, 200],'learning_rate': [0.01, 0.1, 0.2, 1.0],'algorithm': [‘SAMME’, ‘SAMME.R’],'base_estimator__max_depth': [1, 3] algorithm': ‘SAMME.R’, ‘base_estimator__max_depth': 1, ‘learning_rate': 0.1, ‘n_estimators': 100

GBM n_estimators': [50, 100, 200],'max_depth': [3, 5, 7],'learning_rate': [0.01, 0.1, 0.2],'subsample': [0.8, 1.0],'min_samples_split': [2, 5, 10],'min_samples_leaf': [1, 2, 4],'max_features': [‘sqrt’, ‘log2’] learning_rate': 0.1, ‘max_depth': 5, ‘max_features': ‘sqrt’, ‘min_samples_leaf': 1, ‘min_samples_split': 2, ‘n_estimators': 50, ‘subsample': 1.0

SVM kernel': [‘rbf’], ‘C': [0.1, 1, 10, 100],'gamma': [‘scale’, ‘auto’, 0.1, 1, 10] C': 10, ‘gamma': ‘auto’, ‘kernel': ‘rbf'

MLP solver': [‘adam’, ‘sgd’],'learning_rate_init': [0.001, 0.01, 0.1],'momentum': [0.9, 0.95] learning_rate_init': 0.1, ‘momentum': 0.95, ‘solver': ‘sgd'

Models	Hyperparameter	Optimal value
RF	n_estimators’: [50, 100, 200], ‘max_depth': [3, 5, 7, None], ‘max_features': [‘sqrt’, ‘log2’, None],'min_samples_split': [2, 5, 10],'min_samples_leaf': [1, 2, 4],'bootstrap': [True, False],'max_samples': [0.8, 1.0]	bootstrap': True, ‘max_depth': None, ‘max_features’: ‘sqrt’, max_samples': 0.8, ‘min_samples_leaf': 1, min_samples_split': 5, n_estimators': 200
DT	max_depth': [3, 5, 7, None],'min_samples_split': [2, 5, 10],'min_samples_leaf': [1, 2, 4],'max_features': [‘sqrt’, ‘log2’, None],'criterion': [‘gini’, ‘entropy’],'ccp_alpha': [0.0, 0.1, 0.2]	ccp_alpha': 0.0, ‘criterion': ‘gini’, ‘max_depth': 5, max_features': ‘sqrt’, min_samples_leaf': 1, min_samples_split': 5
XGBoost	n_estimators': [50, 100, 200],'max_depth': [3, 5, 7],'learning_rate': [0.01, 0.1, 0.2],'subsample': [0.8, 1.0],'colsample_bytree': [0.8, 1.0]	colsample_bytree': 0.8,'learning_rate': 0.01,'max_depth': 7,'n_estimators': 200,'subsample':0.8
KNN	n_neighbors': [3, 5, 7, 9, 11],'weights': [‘uniform’, ‘distance’],'metric': [‘euclidean’, ‘manhattan’, ‘minkowski’],'p': [1, 2]	metric': ‘manhattan’, ‘n_neighbors': 9, ‘p': 1, ‘weights': ‘uniform'
ET	n_estimators': [50, 100, 200],'max_depth': [3, 5, 7, None],'min_samples_split': [2, 5, 10],'min_samples_leaf': [1, 2, 4],'max_features': [‘sqrt’, ‘log2’, None],'bootstrap': [True, False],'criterion': [‘gini’, ‘entropy’]	bootstrap': False, ‘criterion': ‘entropy’, ‘max_depth': None, ‘max_features': None, ‘min_samples_leaf': 1, ‘min_samples_split': 10, ‘n_estimators': 50
LightGBM	‘num_leaves': [31, 50, 100],'learning_rate': [0.01, 0.05, 0.1],'n_estimators': [100, 200, 500],'min_child_samples': [20, 50, 100],'subsample': [0.6, 0.8, 1.0]	‘num_leaves': 31,'learning_rate': 0.01,'n_estimators': 100,'min_child_samples': 20,'subsample': 0.6
CatBoost	n_estimators': [50, 100, 200],'max_depth': [3, 5, 7], ‘learning_rate': [0.01, 0.1, 0.2],'subsample': [0.8, 1.0],'colsample_bylevel': [0.8, 1.0],'border_count': [32, 64, 128],'bagging_temperature': [0, 1, 2],'l2_leaf_reg': [0, 1, 3]	bagging_temperature': 0, ‘border_count': 128, ‘colsample_bylevel': 1.0, ‘l2_leaf_reg': 0, ‘learning_rate': 0.2, ‘max_depth': 3, ‘n_estimators': 200, ‘subsample': 1.0
AdaBoost	n_estimators': [50, 100, 200],'learning_rate': [0.01, 0.1, 0.2, 1.0],'algorithm': [‘SAMME’, ‘SAMME.R’],'base_estimator__max_depth': [1, 3]	algorithm': ‘SAMME.R’, ‘base_estimator__max_depth': 1, ‘learning_rate': 0.1, ‘n_estimators': 100
GBM	n_estimators': [50, 100, 200],'max_depth': [3, 5, 7],'learning_rate': [0.01, 0.1, 0.2],'subsample': [0.8, 1.0],'min_samples_split': [2, 5, 10],'min_samples_leaf': [1, 2, 4],'max_features': [‘sqrt’, ‘log2’]	learning_rate': 0.1, ‘max_depth': 5, ‘max_features': ‘sqrt’, ‘min_samples_leaf': 1, ‘min_samples_split': 2, ‘n_estimators': 50, ‘subsample': 1.0
SVM	kernel': [‘rbf’], ‘C': [0.1, 1, 10, 100],'gamma': [‘scale’, ‘auto’, 0.1, 1, 10]	C': 10, ‘gamma': ‘auto’, ‘kernel': ‘rbf'
MLP	solver': [‘adam’, ‘sgd’],'learning_rate_init': [0.001, 0.01, 0.1],'momentum': [0.9, 0.95]	learning_rate_init': 0.1, ‘momentum': 0.95, ‘solver': ‘sgd'

3.3.2 Model performance comparison

As can be seen from Table 3, the CAT model shows the best performance in experiments with a 20% test set share, with a classification correctness and recall of 81.19% and 80.00%, respectively, and its F1 score (0.7912) and NPV (0.8364) are in the top rankings, which indicates a strong discriminative ability for positive and negative class samples. The model demonstrated a well-balanced performance in the classification task and was able to effectively identify the risk of diabetic retinopathy (DR).The LGBM model had a high specificity of 83.93% (NPV = 0.8174), which was outstanding in controlling the false positives and is suitable for scenarios sensitive to the cost of misdiagnosis.The XGB and the ET model performed close to each other in terms of accuracy and sensitivity, but ET's F1 score (0.7692) is slightly lower than that of XGB (0.7760), reflecting that XGB may be better adapted under category-imbalanced data.The recall/sensitivity (≤71.11%) of DT and KNN is significantly low, which may lead to more DR cases being missed and higher clinical risk. The multilayer perceptron (MLP) and support vector machine (SVM), although the accuracy was up to standard, the specificity (≤79.46%) was insufficient, which may increase the misdiagnosis rate of non-DR patients.

Table 3.
Model performance results with different test set occupancies.

Test set Models Accuracy Recall F1 Score Sensitivity Specificity PPV NPV

20% RF 0.7871 0.7667 0.7624 0.7667 0.8036 0.7582 0.8108

DT 0.7525 0.7111 0.7191 0.7111 0.7857 0.7273 0.7719

XGB 0.7970 0.7889 0.7760 0.7889 0.8036 0.7634 0.8257

KNN 0.7673 0.7000 0.7283 0.7000 0.8214 0.7590 0.7731

ET 0.7921 0.7778 0.7692 0.7778 0.8036 0.7609 0.8182

LGBM 0.8069 0.7667 0.7797 0.7667 0.8393 0.7931 0.8174

CAT 0.8119 0.8000 0.7912 0.8000 0.8214 0.7826 0.8364

ADA 0.7970 0.7444 0.7657 0.7444 0.8393 0.7882 0.8034

GBM 0.7822 0.7222 0.7471 0.7222 0.8304 0.7738 0.7881

SVM 0.7970 0.8000 0.7784 0.8000 0.7946 0.7579 0.8318

MLP 0.7871 0.7778 0.7650 0.7778 0.7946 0.7527 0.8165

30% RF 0.8079 0.7836 0.7836 0.7836 0.8274 0.7836 0.8274

DT 0.7748 0.6642 0.7236 0.6642 0.8631 0.7946 0.7632

XGB 0.8013 0.7537 0.7710 0.7537 0.8393 0.7891 0.8103

KNN 0.8013 0.7836 0.7778 0.7836 0.8155 0.7721 0.8253

ET 0.8245 0.8209 0.8059 0.8209 0.8274 0.7914 0.8528

LGBM 0.8046 0.7612 0.7757 0.7612 0.8393 0.7907 0.8150

CAT 0.8079 0.7761 0.7820 0.7761 0.8333 0.7879 0.8235

ADA 0.7715 0.7090 0.7336 0.7090 0.8214 0.7600 0.7797

GBM 0.8179 0.7761 0.7909 0.7761 0.8512 0.8062 0.8266

SVM 0.8113 0.8284 0.7957 0.8284 0.7976 0.7655 0.8535

MLP 0.7815 0.8134 0.7676 0.8134 0.7560 0.7267 0.8355

Test set	Models	Accuracy	Recall	F1 Score	Sensitivity	Specificity	PPV	NPV
20%	RF	0.7871	0.7667	0.7624	0.7667	0.8036	0.7582	0.8108
DT	0.7525	0.7111	0.7191	0.7111	0.7857	0.7273	0.7719
XGB	0.7970	0.7889	0.7760	0.7889	0.8036	0.7634	0.8257
KNN	0.7673	0.7000	0.7283	0.7000	0.8214	0.7590	0.7731
ET	0.7921	0.7778	0.7692	0.7778	0.8036	0.7609	0.8182
LGBM	0.8069	0.7667	0.7797	0.7667	0.8393	0.7931	0.8174
CAT	0.8119	0.8000	0.7912	0.8000	0.8214	0.7826	0.8364
ADA	0.7970	0.7444	0.7657	0.7444	0.8393	0.7882	0.8034
GBM	0.7822	0.7222	0.7471	0.7222	0.8304	0.7738	0.7881
SVM	0.7970	0.8000	0.7784	0.8000	0.7946	0.7579	0.8318
MLP	0.7871	0.7778	0.7650	0.7778	0.7946	0.7527	0.8165
30%	RF	0.8079	0.7836	0.7836	0.7836	0.8274	0.7836	0.8274
DT	0.7748	0.6642	0.7236	0.6642	0.8631	0.7946	0.7632
XGB	0.8013	0.7537	0.7710	0.7537	0.8393	0.7891	0.8103
KNN	0.8013	0.7836	0.7778	0.7836	0.8155	0.7721	0.8253
ET	0.8245	0.8209	0.8059	0.8209	0.8274	0.7914	0.8528
LGBM	0.8046	0.7612	0.7757	0.7612	0.8393	0.7907	0.8150
CAT	0.8079	0.7761	0.7820	0.7761	0.8333	0.7879	0.8235
ADA	0.7715	0.7090	0.7336	0.7090	0.8214	0.7600	0.7797
GBM	0.8179	0.7761	0.7909	0.7761	0.8512	0.8062	0.8266
SVM	0.8113	0.8284	0.7957	0.8284	0.7976	0.7655	0.8535
MLP	0.7815	0.8134	0.7676	0.8134	0.7560	0.7267	0.8355

When the test set share was increased to 30%, the ET model performed outstandingly on the 30% test set, achieving the best balance of 82.45% accuracy and 82.09% sensitivity/specificity, and its NPV (0.8528) indicated excellent ability to correctly exclude non-DR patients.The ET model still maintained high classification accuracy and stability on the larger test set, indicating its ability to adaptability to data distribution.The GBM model has a high specificity of 85.12% (NPV=0.8266) and performs optimally in controlling false positives, making it suitable for misdiagnosis-sensitive clinical scenarios.The SVM model has the highest recall/sensitivity (82.84%) and the strongest ability to detect DR cases, but has a low specificity (79.76%), which needs to be weighed against the risk of missed and misdiagnosed cases. In addition, RF, GBM and SVM models also showed better performance, with SVM achieving 82.84% in recall, demonstrating a stronger ability to identify positive class samples (high-risk patients).DT still had the lowest recall (66.42%) and sensitivity (66.42%), with a significant risk of clinical underdiagnosis.The specificity (≤75.60%) and the ADA (≤75.60%) of the MLP and the PPV (≤76.00%) were low, which may increase the misdiagnosis rate in non-DR patients.

The generalization performance of the models can be observed by comparing the results of the 20% and 30% test sets as shown in Table 3: the ET model performs robustly on both test sets, with fluctuations of accuracy and sensitivity less than 2%, indicating that its predictive ability is less affected by the sample partitioning.The CAT model performs well on the 20% test set, but the sensitivity decreases by 2.39% on the 30% test set, which may be due to the increased amount of data exposing the risk of overfitting. The XGB and LGBM models have less fluctuation in performance, showing their robustness under large samples.The SVM model has an increase in recall by 2.84% but a decrease in specificity by 2.7% at 30% test set, reflecting that its discriminative ability for positive and negative class samples is sensitive to the change in data volume.

3.3.3 Exploring fusion machine learning algorithms

VotingClassifier is an integrated learning approach that combines several different machine learning models to make decisions about the classification task by voting, thus improving the predictive performance and robustness of the overall model.VotingClassifier provides two main voting strategies hard voting and soft voting, hard voting is where each submodel predicts for each sample a category label, and the final classification result depends on the category with the most votes for the predicted category in the sub-model. In this paper, we implement a soft voting mechanism (Voting_soft) using dynamically weighted integrated learning based on model confidence. The integration weights are dynamically assigned by quantifying the uncertainty (entropy) of the prediction probability of each base model: firstly, the Shannon entropy is calculated for the prediction probability of each base model on the validation set (a lower entropy value indicates that the prediction probability is more concentrated, and the higher the confidence is), and the inverse of the entropy is used as the initial weights in order to avoid that the low confidence model dominates the integration results, and finally the probability weights are obtained through the normalization process. These weights are used to construct soft-voting classifiers to make the prediction results of high-confidence models play a greater role in integration decisions. This approach breaks through the limitations of traditional static weighting or simple majority voting, and enables a smarter model fusion strategy by dynamically capturing the differences in model reliability on specific tasks. The base model here consists of RandomForest, XGBoost, CatBoost, and ExtraTrees in subsection 3.3.2.

The ROC curves are plotted against the PR curves for a 20% test set share in Figure 2 A and Figure 2B. The ROC curves are plotted against the PR curves for a 30% test set share in Figure 2C and Figure 2D. The corresponding confidence intervals for the corresponding values are given. The 95% confidence intervals corresponding to the corresponding values are also given. It can be seen that the AUC values of the fusion machine learning model constructed by the soft voting mechanism proposed in this paper are optimal in the two datasets with different occupancy ratios, which are 0.8911 and 0.9040, respectively.In the comparison of the area AR values under the PR curves, Voting_soft (0.8558) is only slightly lower than the Random Forest (0.8580) under the 20% testing set occupancy ratio, and Voting_soft (0.8558) is only slightly lower than the Random Forest (0.8580) under the 20% testing set occupancy ratio. 20% test set share Voting_soft (0.8580) is only slightly lower than XGBoost (0.8659), which fully demonstrates the advantage of Voting_soft model in diabetic retinopathy risk prediction.

Figure 2.

Graph of ROC curve and PR curve.

3.4 Assessment of DR impact factors

In order to deeply analyze the contribution of each risk factor to the prediction of diabetic retinopathy (DR), the study used the optimal XGBoost model combined with the SHAP interpretation framework to carry out the analysis. Figure 3 systematically demonstrates the results of the model interpretability analysis based on SHAP values: the framework reveals the association mechanism between clinical indicators and the development of DR through the quantification of feature attribution, in which the positive and negative directions of the SHAP values directly correspond to the risk increasing and decreasing effects, the vertical distribution of the sample points reflects the feature influence weight, and the color gradient characterizes the size of the indicator values. The global feature importance ranking in Figure 3A indicates that the indicators of AGE, TG, SEX and HDLC constitute the core decision factors for DR prediction. The SHAP summary plot in Figure 3B further validates the above findings by visualizing the marginal contribution of each risk factor to model prediction through scatter distribution, and specifically points out that TG metrics present a significant positive driving effect - samples with high TG levels (red dot set) are concentrated in high-risk regions, while samples with low TG values (blue dots set) correspond more to risk suppression effects.

Figure 3.

SHAP analysis diagram.

In response to the unique SHAP value calculation characteristics of the XGBoost model, the study innovatively disassembles the analysis by real labels: Figure 3C and Figure 3D construct the box plots of DR negative (class 0) and positive (class 1), respectively, to realize the differentiated interpretation of the prediction results of the two classes. The force-directed plot in Figure 3E visualizes the superposition process of the model baseline value and the sample characteristic contribution, and it should be noted in particular that the f(x) value in the plot represents the log odds ratio before Sigmoid transformation, rather than the direct probability output. Figure 3F also shows the impact of each feature on the predictive performance of the model, reflecting the magnitude of the SHAP value for each feature. The blue line in the figure indicates the trend of the AUC with respect to the features, and the yellow bar highlights the features that have the greatest impact on the predictive performance. With this visualization, it can be seen that the AUC of the XGB model stabilizes to more than 90% after incorporating the 7th most important factor, indicating that AGE, TG, SEX, HDLC, PR, BUN, and VLDLC are important risk factors for DR, especially the first two. This progressive feature utility analysis not only verified the strong association between metabolic indicators and DR, but also provided a quantifiable risk factor prioritization basis for clinical intervention strategies.

4 Discussion

This study compared multiple machine learning algorithms using patient data and identified a set of predictive factors for diabetic retinopathy (DR). Model performance was evaluated using accuracy, precision, recall, F1-score, ROC curves, and PR curves. Many prior studies have used machine learning to predict DR. Table 4 compares this research with existing literature. Hong DD²⁸ recruited 1739 type 2 diabetes patients and constructed DR prediction models using five algorithms (e.g., logistic regression, decision tree). Results showed increased DR risk at PR values of 57 and 72 beats/min, and elevated risk with rising BUN levels (≥5 mmol/L), consistent with our SHAP analysis. Zhao et al.²⁹ conducted a retrospective study on 7943 Chinese patients, using RF-based feature selection to identify predictors. Among RF, XGBoost, LR, SVM, and KNN models, XGBoost performed best (accuracy: 80.30%, AUC: 0.889). Tsao et al.³⁰ proposed an SVM-based model (accuracy: 79.0%, AUC: 0.839), outperforming LR, ANN, and DT. Ogunyemi et al.³¹ analyzed 24 predictors in 513 U.S. patients, using LASSO for feature selection. The RUSBoost model achieved 73.5% accuracy. Oh et al.³² applied ridge regression, LASSO, and elastic net to predict DR in 490 subjects. Hosseini et al.³³ included 3734 Iranian patients with 11 predictors, using LR-based diagnosis without feature selection. In our patient dataset, the LightGBM model outperformed others. SHAP analysis identified TG, HDLC, and BUN as primary predictors, consistent with findings in.²⁸^33–35 Additionally, AGE and SEX were identified as key factors, with DR risk increasing significantly with age, aligning with literature.⁶^36–38

Table 4.
Comparison of the proposed study with the existing studies in literature.

Author's Year Data size Feature selection method Accuracy(%) SHAP

Hong DD²⁸ 2024 1739 Random Forest 79.70 YES

Zhao et al.²⁹ 2022 7943 Random Forest 80.30 NO

Tsao et al.³⁰ 2018 536 - 79.50 NO

Ogunyemi et al.³¹ 2015 513 LASSO 73.05 NO

Oh et al.³² 2013 490 - 75.20 NO

Hosseini et al.³³ 2009 3734 - 70.40 NO

In this paper 1006 Random Forest 80.46 YES

Author's	Year	Data size	Feature selection method	Accuracy(%)	SHAP
Hong DD²⁸	2024	1739	Random Forest	79.70	YES
Zhao et al.²⁹	2022	7943	Random Forest	80.30	NO
Tsao et al.³⁰	2018	536	-	79.50	NO
Ogunyemi et al.³¹	2015	513	LASSO	73.05	NO
Oh et al.³²	2013	490	-	75.20	NO
Hosseini et al.³³	2009	3734	-	70.40	NO
In this paper		1006	Random Forest	80.46	YES

5 Conclusion and outlook

This study delves into the application of interpretable machine learning algorithms in the diagnostic prediction of DR and significant research results have been achieved. By comparing the performance of different machine learning algorithms, we found that models such as LightGBM, XGBoost and Random Forest perform well in DR risk prediction and significantly outperform traditional diagnostic methods. In particular, the LightGBM model shows the best performance when dealing with complex medical datasets, with high AUC values, highlighting its advantages in improving prediction accuracy. In terms of model interpretability, the SHAP method provides an effective tool to explain and understand the decision-making process of machine learning models. By analyzing the SHAP values, we found that characteristics such as age, TG, SEX, and HDLC played a key role in predicting the prevalence of DR. This finding not only enhances the clinical application value of the model, but also helps physicians to better understand the model decisions and improve clinical trust and acceptance.

Although this paper has achieved certain results in DR diagnosis, the research is based on a single data source. Future research will integrate multimodal models of computer vision and clinical structured data to significantly improve the accuracy and comprehensiveness of DR risk prediction. Moreover, although the model performs well on a specific dataset, its ability to generalize still needs to be improved. Future research will attempt to adopt techniques such as transfer learning, semi-supervised learning, and deep learning to improve the applicability and accuracy of models in different populations and environments. At the same time, strengthening the research on model interpretability is also an important direction for future development, which will help support clinical decision-making and the development of personalized medicine.

Footnotes

Acknowledgments

Data availability statement: The data are available from the corresponding author upon request.

Ethical statement

None to report.

Author contributions

Yifeng Dou and Jiantao Liu conceived the study. Yifeng Dou participated in the study design, data analysis and statistics. Yifeng Dou and Jiantao Liu helped draft the manuscript. All authors read and approved the final manuscript.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the the In-hospital project of Tianjin Medical University Baodi Hospital, (grant number BDYYQN01).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

The data are available from the corresponding author upon request.

ORCID iDs

Yifeng Dou

Jiantao Liu

References

Yin

. Diabetic retinopathy: looking forward to 2030. Front Endocrinol (Lausanne). 2023; 13: 1077669.

Lin

Hsih

Lin

, et al. Update in the epidemiology, risk factors, screening, and treatment of diabetic retinopathy. J Diabetes Investig. 2021; 12: 1322–1325.

Ling

YihChung

Chak

MYY

, et al. Global prevalence of diabetic retinopathy and projection of burden through 2045: systematic review and meta-analysis. Ophthalmology 2021; 128: 1580–1591.

Shi

. Progress of machine learning in clinical research. China Digit Med 2025; 20: 1–10.

Sun

Liu

, et al. Research progress of deep learning in the field of diabetic retinopathy classification. Comput Eng Applic. 2024; 60: 16–30.

Wang

. Research on the risk prediction model of diabetic retinopathy occurrence based on machine learning. Qufu Normal University, 2024.

. Risk factor analysis and prediction model construction of diabetic retinopathy. Jilin University, 2024.

Kim

Park

Son

, et al. Development and validation of a machine learning algorithm for predicting diabetes retinopathy in patients with type 2 diabetes: algorithm development study. JMIR Med Inform. 2025; 13: e58107.

Luong

Cheung

McMurtry

, et al. Comparison of machine learning models to a novel score in the identification of patients at low risk for diabetic retinopathy. Ophthalmol Sci 2025; 5: 100592–100592.

10.

Zhao

Jin

Xiao

, et al. Application of a multimodal model integrating computer vision and structured data in diabetic retinopathy referral. China Digit Med 2024; 19: 29–35.

11.

Zhang

, et al. The application and clinical translation of the self-evolving machine learning methods in predicting diabetic retinopathy and visualizing clinical transformation. Front Endocrinol (Lausanne). 2024; 15: 1429974–1429974.

12.

Pang

Luo

Zhang

, et al. Multi-Omics integration with machine learning identified early diabetic retinopathy, diabetic macula edema and anti-VEGF treatment response. Transl Vis Sci Technol. 2024; 13: 23.

13.

Lukashevich

. Diabetic retinopathy fundus image classification using ensemble methods. Pattern Recognit Image Anal. 2024; 34: 331–339.

14.

Abdalla

IMM

Mohanraj

. Revolutionizing diabetic retinopathy screening and management: the role of artificial intelligence and machine learning. World J Clin Cases. 2025; 13: 101306.

15.

Breiman

. Random forests. Mach Learn. 2001; 45: 5–32.

16.

Mitchell

. Machine learning . 1997.

17.

Chen

Benesty

. xgboost: Extreme Gradient Boosting. 2016.

18.

Kramer

. K-Nearest neighbors. Berlin Heidelberg: Springer, 2013.

19.

Geurts

Ernst

Wehenkel

. Extremely randomized trees. Mach Learn. 2006; 63: 3–42.

20.

Meng

. LightGBM: A highly efficient gradient boosting decision tree. In Neural information processing systems. Curran Associates Inc., 2017.

21.

Wang

Huang

Ren

, et al. In-process belt-image-based material removal rate monitoring for abrasive belt grinding using CatBoost algorithm. Int J Adv Manuf Technol 2022; 123: 2575–2591.

22.

Morra

Apostolova

, et al. Comparison of AdaBoost and support vector machines for detecting Alzheimer's disease through automated hippocampal segmentation. IEEE Trans Med Imaging. 2010; 29: 30–43.

23.

Rokach

. Gradient boosting machines. Series in Machine Perception and Artificial Intelligence, 2019.

24.

Liu

Zhu

Dai

, et al. Multi-task support vector machine classifier with generalized huber loss. J Classif. 2025; 42: 221–252.

25.

Sivakumar

Desai

. Image restoration using a multilayer perceptron with a multilevel sigmoidal function. IEEE Trans Signal Process. 1993; 41: 2018–2022.

26.

Zhao

Wang

, et al. A machine-learning-derived online prediction model for depression risk in COPD patients: a retrospective cohort study from CHARLS. J Affect Disord. 2025; 377: 284–293.

27.

Lundberg

Gabriel

Hugh

, et al. From local explanations to global understanding with explainable AI for trees. Nat Mach Intell 2020; 2: 56–67.

28.

Hong

. Investigation on the prevalence of diabetic retinopathy and construction of clinical prediction model in plain wind-sand and loess hilly areas of Gansu Province. Gansu University of Chinese Medicine, 2024.

29.

Zhao

, et al. Using machine learning techniques to develop risk prediction models for the risk of incident diabetic retinopathy among patients with type 2 diabetes mellitus: a cohort study. Front Endocrinol 2022; 13: 876559.

30.

Tsao

Chan

. Predicting diabetic retinopathy and identifying interpretable biomedical features using machine learning algorithms. BMC Bioinf 2018; 19: 111–121.

31.

Ogunyemi

Kermah

. Machine learning approaches for detecting diabetic retinopathy from clinical and public health records. AMIA Annu Symp Proc 2015; 2015: 983–990.

32.

Yoo

Park

. Diabetic retinopathy risk prediction for fundus examination using sparse learning: a cross-sectional study. BMC Med Inf Decis Making 2013; 13: 1–4.

33.

Hosseini

Maracy

Amini

, et al. A risk score development for diabetic retinopathy screening in Isfahan-Iran. J Res Med Sci Off J Isfahan Univ Med Sci 2009; 14: 105.

34.

Islam

Rahman

Rabby

, et al. Predicting the risk of diabetic retinopathy using explainable machine learning algorithms. Diabetes Metab Syndr Clin Res Rev 2023; 17: 11.

35.

Guan

. Analysis of risk factors and diagnostic model for diabetic nephropathy in patients with diabetic retinopathy. Anhui Medical University, 2024.

36.

Goldstein

Ding

Carasquillo

, et al. Prediction of proliferative diabetic retinopathy using machine learning in Latino and non-Hispanic black cohorts with routine blood and urine testing. Ophthalm Physiol Opt J Br College Ophthalm Optic (Optometrists) 2024; 45: 1549–1567.

37.

Yang

. Construction of clinical prediction model for retinopathy in Type 2 diabetes. Lanzhou University, 2023.

38.

Tebeje

Yenit

Nigatu

, et al. Prediction of diabetic retinopathy among type 2 diabetic patients in University of Gondar Comprehensive Specialized Hospital, 2006–2021: a prognostic model. Int J Med Inf. 2024; 190: 105536.