Early Prediction of Gestational Diabetes Using Deep Learning and Electrocardiography: Results from the MERGD Study

1. Introduction

Gestational diabetes (GD) continues to be an important antenatal complication that can adversely affect maternal and fetal outcomes. The International Diabetes Federation (IDF) Diabetes Atlas estimated that the global prevalence of GD is as high as 14% of pregnancies.[1] However, there exists wide variation in the prevalence estimates for GD across continents, across countries and states within countries.[1] A primary reason for this variation is a lack of unified definition of GD. For example, the American College of Obstetrics and Gynecology (ACOG) recommends a two-step procedure for GD diagnosis[2] – the first step is a glucose challenge test (GCT) with 50g glucose, and the second step is an oral glucose tolerance test (OGTT). The GCT involves a measurement of blood glucose 1 hour after ingestion of the glucose load and patients with ≥130 mg/dl (≥7.2 mmol/L) of blood glucose are recalled at 24-28 weeks for a 3-hour OGTT. This method benefits from the combination of a sensitive (GCT) and a specific (OGTT). However, in low- to middle-income countries, it is pragmatically challenging to conduct a two-step procedure for GD diagnosis.[3,4,5] Therefore, different countries adopt different practices and definitions of GD diagnosis.[6] Unification of such diverse approaches to arrive at standardized, comparable estimates of GD prevalence is therefore needed. While OGTT remains the confirmatory test, the step of glucose challenge test is less suited to low-resource settings.

Considering this, several attempts have been made to develop and validate early risk assessment strategies for gestational diabetes. For example, the Monash Early Pregnancy Risk Score,[7] uses clinical variables such as age, height, weight, ethnicity, family history of diabetes and past history of GD to predict subsequent risk of GD during the first trimester of pregnancy. Similarly, the Tianjin Women and Children's Health Center Model, validated in large Chinese cohorts additionally uses information on physical activity, passive smoking and weight gain.[8] The accuracy of these scoring systems was 70.2% and 71.2%, respectively. The pre-pregnancy biomarker risk score developed by Badon et al[9] used biomarkers such as sex hormone-binding globulin, adiponectin, fasting glucose and homeostatic model for assessment of insulin resistance to achieve an accuracy of 73%. The recently developed gestational diabetes risk prediction model[10] achieved an improved accuracy of 78% by additionally including glycated hemoglobin (HbA1c).In India, a recent study that combined clinical and biochemical predictors achieved an accuracy of 92% when used during the second trimester.[11] Together these studies indicate that improved and more accurate models are needed to predict the risk of GDM early during pregnancy.

There is now burgeoning evidence to support the hypothesis that electrocardiographic (ECG) changes accompany glucose level variations and that ECG patterns are characteristic of prediabetes and diabetes.[5,12,13] Furthermore, ECG changes have been described in women who developed gestational diabetes.[14] We, therefore, hypothesized that the potential diagnosis of GD may be reflected in subtle ECG changes early during pregnancy. If this hypothesis is supported by data, then ECG can be potentially used as a first screen at the first possible healthcare contact during pregnancy. However, such subtle and non-linear changes in ECG are unlikely to be picked up using clinical gestalt or traditional association statistics. Therefore, we conducted this study where we combined the information captured by the ECG with novel feature extraction procedures and the power of deep learning to discriminate identify women who, on follow-up, progressed to GD.

2. Materials and Methods

2.1. Study Participants

Data for this study came as a part of the Markers of Early Risk-stratification of Gestational Diabetes (MERGD) study (registered with the Clinical Trials Registry—India, CTRI/2018/05/013946).The study protocol has been described in detail previously.[15] Briefly, the MERGD study was conducted on eligible and consenting pregnant women who reported to the Daga Memorial Women’s Hospital, Nagpur – a secondary care hospital specializing in obstetric care. Inclusion criteria were consecutive, newly registered pregnant women at the Study Center, with a gestational age at first contact < 20 weeks, no history of type 2 diabetes, and who had provided written, informed consent. Enrollment took place between 21 May 2018 and 11 August 2018 and participants were followed up for outcomes assessment till the end of delivery. Last date of follow-up was 22 February 2019. This study was approved by the Ethics Research Committee of the Daga Memorial Women’s Hospital, Nagpur, India on 10 May 2018. A total of 1,040 women were enrolled in the study.

2.2. Ground Truth

Gestational diabetes was defined as a compound phenotype as follows: The study used a two-step procedure recommended by the American College of Obstetrics and Gynecology (ACOG) in which the first step was a 50 g glucose load, 1 h glucose challenge test (GCT), where a cutoff post-load plasma glucose concentration of 130 mg/dL (7.2 mmol/L) was used to decide the need for an oral glucose tolerance test (OGTT). The second step was carried out in women with an abnormal GCT value and included a 100 g glucose load OGTT using a 100 g glucose load for a 3 h OGTT. Abnormal glucose values were investigated at the time of glucose load (fasting) and then hourly after glucose load. An abnormal value was defined using both the Carpenter–Coustan (C&C) criteria[16] and the National Diabetes Data Group (NDDG) criteria.[17] The C&C criteria used were: fasting—≥95 mg/dL (5.28 mmol/L); 1 h—≥180 mg/dL (10.0 mmol/L); 2 h—≥155 mg/dL; and 3 h—≥140 mg/dL (7.78 mmol/L). The NDDG criteria were: fasting—≥105 mg/dL, one hour—≥190 mg/dL, two hour—≥165 mg/dL and three hour—≥145 mg/dL. Presence of one or more abnormal values was defined as GD. In addition, an HbA1c of ≥6.5% was also considered as GD.

2.3. ECG Recording and Preprocessing

Study participants underwent the standard, 12-lead ECG which was recorded using a digital ECG device (Cardiart 6208, BPL Medical Technologies, Mumbai, Maharashtra, India). All recordings were for a duration of 10 s, with ≥12 bit analogue to digital converters, 0.05–150 Hz bandwidth and a sampling frequency of 1000 Hz. ECG recordings were saved as dicom files and deidentified before analytical use.

Considering the sampling frequency and duration, each ECG resulted in an array of 10,000x12 readings. We then used the Neurokit 2[18] Python library for initial preprocessing of the ECG. Each lead was separately preprocessed using the ecg.clean() analytical pipeline. Under this, the lead-specific ECG corrected for baseline wander and then followed the Pan & Tompkins algorithm.[19] This included bandpass filtering (to isolate 5-40 Hz frequency, QRS detection using differentiation, sliding widow integration and adaptive thresholding. Further, we used the Zhao method[20] to determine the ECG quality.

The cleaned, lead-specific ECG signals were then stacked back to into a two-dimensional array of 10,000x12 array. The full dataset thus represented an array of size 1,036x10,000x12 as shown in Figure 1. We then used the SciPy library[21] to resample the signal at 100 Hz and therefore the final data array had a size of 1,036x1000x12. This array was saved as a Numpy[22] binary file and used for subsequent analyses.

2.4. Feature Extraction

Our first step was to extract features from the ECG signal. We extracted features from three sources – using an autoencoder, collating basic time domain characteristics reported by the recording device and extracting features of heart rate variability (HRV) extracted by the Neurokit 2 library.

2.4.1. Autoencoder-Based Feature Extraction

The motivation for this step was to represent the entire 10s ECG signal for each lead in as few but informative number of features as possible. For this, we used an autoencoder framework which comprised of the encoder and decoder components as shown in Figure 1. The encoder used a series of compressing blocks each of which contained a convolution neural network (CNN) with prespecified number of 3x3 kernels, a batch normalization layer and a layer that applied the leaky rectified unit (leaky relu) activation function. The number of CNN layers started with 1000 (the input dimension) and proceeded to 512, 256, 128 and 64. This last layer of 64 kernels represented the bottleneck layer and the decoder component of the autoencoder used this layer as input. The decoder proceeded in a mirror image fashion and recreated the 1000 points that the encoder started out with. The loss function for the autoencoder was the mean squared error (MSE) that compared the input of the encoder function to the output of the decoder layer. The autoencoder was trained using tensorflow/keras platform.[23,24] Since we aimed to characterize a normal ECG signal, we used data only from women who did not develop GD for the autoencoder training. The autoencoder model was trained on a 2/3 training set and validated on the remaining 1/3 data. The best autoencoder model that provided the minimum MAE was used for feature extraction. To extract features, we predicted each lead of each 12-lead ECG into a 1x64 feature map drawn from the output of the encoder component. Thus, for each ECG we had 12 such feature maps which were finally horizontally stacked to create a feature vector of dimension 1x768 and contained a total of 768 features extracted from all 12 leads. Together, the feature array was of size 1,036x768 representing 768 features extracted from each ECG.

2.4.2. Basic Features Reported by the ECG Device

The ECG recording machine reported basic time domain features for each recorded ECG at the time of recording. These features (N_features = 10) were heart rate, duration of the P wave, PR interval, duration of the QRS complex, QT_corrected interval, mean axis for the P, QRS and T waves and two voltages – the R peak in lead V5 and the S depth in lead V1. From these we created an array of size 1,036x10 as the second set of features.

2.4.3. Features Related to Heart Rate Variability

The Neurokit 2 library has an in-built function (ecg.hrv()) to estimate several features related to heart rate variability from a given ECG signal. We used the signal from lead II of each ECG to estimate these features (N_features = 40) and generated a 1,036x40 array of HRV-based features. Lastly, we horizontally stacked the features arrays from each of three abovementioned sources to construct a full feature map of size 1,036x818 (represented as a black rectangle in Figure 1).

2.5. Generation of Synthetic Data for Classification Model Training

To develop and train a classification model, we aimed to create a synthetic dataset based off the actual feature map (1,036x818) described above. For this, we created 15 copies of the original dataset each of which was augmented by adding Gaussian noise with different parameters. Five values of mean (-0.2, -0.1, 0, 0.1, and 0.2) and three values of standard deviation (0.1, 0.2, and 0.3) yielded these 15 combinations. Each copy had the same number of women who developed GD. In the next step, we used the synthetic minority oversampling technique (SMOTE)[25] to balance each copy and then vertically stacked these datasets to generate an array with 27,690 rows and 818 columns. This array represented the synthetic dataset which was then minmax standardized and used for further analyses. It should be noted that not a single instance of the 27,690 rows in this dataset directly represented any instance in the original data of 1,036x818 features. Thus, the model was trained on a completely synthetic dataset and then validated on the original, unseen data.

2.6. Random Forest Model Training

We used the Scikit-Learn library[26] to implement, train, cross-validate and then internally validate a random forest (RF) model for classification. We first split the synthetic dataset into a training and test set (95% and 5%, respectively. The training set was used for model development. Best performing model was chosen from a combination of hyperparameter grid search and the 10 folds based on the training set. Following hyperparameter space was used to generate different combinations in the grid search: maximum tree depth (2,3,4), maximum number of features (2,3,4), minimum number of samples per leaf (1,2,3), minimum samples per split (2,3,4) and number of trees (75,100). Together these hyperparameter choices and the 10 folds meant that a total of 1,620 model fits were assessed to choose the best performing model. The loss function used is optimization process was logloss for binary classification.

2.7. Validation

We validated the RF model in three steps on three different datasets. First, we used 10-fold cross validation during model training to ensure that the model performed consistently within the training set. Second, we chose the best model to assess its classification performance in the test set. Third, we evaluated the classification performance of the best model in the original dataset of 1,036x818 features in two ways. We first used receiver operating characteristic (ROC) curve-based statistics to quantitatively measure the predictive performance of the model and then we compared this performance with that of clinical predictors only. Further, we compared the predictive performance of the RF model with that of well-established, externally validated risk scoring systems used to predict the risk of GD. A recent, comprehensive meta-analysis[27] recommended four risk scoring systems by Naylor et al,[28] Teede et al,[7] van Leeuwen et al[29] and Nanda et al[30] as early screens for the risk of GD. We compared the performance of our RF model with these four risk scoring systems as well.

2.8. Statistical Analyses

Descriptive statistics included mean and standard deviation for continuous variables and number and percentages for categorical variables. Difference between women with and without GD was tested using Student’s t test for continuous variables, and χ² test for categorical variables. In case of cell counts less than 5 for categorical variables, Fisher’s exact test was used in place of the χ² test. Classification performance was assessed by plotting the ROC curve and estimating the area under the ROC curve (AUC). Comparison of AUC of two curves was done using the De Long and De Long test. Best threshold to evaluate a test performance was chosen as the point on the ROC curve nearest to the upper-left corner. Dichotomous classification performance was assessed by estimating sensitivity (true positives/total positives), specificity (true negatives/total negatives), precision (true positives/predicted positives), accuracy (proportion of correctly classified individuals) and F1 score (harmonic mean of sensitivity and precision). We also estimated the likelihood ratio of positive (LR+) and negative (LR-) prediction of GD. Importance of individual features was assessed using mean decrease in impurity (MDI) metric[31] and informative number of features was estimated using the kneedle algorithm.[32] Statistical analyses were conducted using dedicated python scripts and using the Stata 18.0 (Stata Corp, College Station, TX) software package. Statistical significance was tested at a global type I error rate of 0.05.

3. Results

3.1. Study Participants

Details of the MERGD cohort have been described previously.[15] In this study, we included a total of 1,036 participants for whom first visit ECGs were available. Table 1 describes and compares the baseline characteristics of participants who did (n = 123, 11.87%) and did not (88.13%) develop GD. Briefly, women who developed GD were, on an average, 1.31 years older and registered 0.65 weeks earlier for antenatal care as compared to the women who did not develop GD. Family income, caste and religion were all associated with the risk of GD. Specifically, of the women who developed GD, 20.83% belonged to families earning ≥ ₹ 200,000 per annum, 44.72% belonged to the open caste and 32.52% were Muslim.

Comparatively, these numbers in the women who did not develop GD were 13.41%, 30.04% and 32.52%, respectively and were statistically significant (p = 0.0291, 0.0010, and 0.0029, respectively). Further, women who developed GD had a significantly higher body mass index (BMI, p = 0.0003), systolic blood pressure (p = 0.0012), diastolic blood pressure (p = 0.0001) and mean arterial pressure (p = 0.0001). Lipid profile showed that women with GD had significantly higher serum triglycerides and very low-density lipoproteins (p = 0.0063 each).

As per the MERGD study protocol, the date of conception was determined using ultrasonography (USG, Figure 2). Using this date as day 0, the average gestational age at the first visit was 16.15 weeks, a glucose challenge test (GCT) was conducted at an average of 24.66 weeks and oral glucose tolerance test (OGTT, if needed based on GCT result) at 26.56 weeks. Thus, there was a lead time of 10.41 weeks between ECG measurement and OGTT-confirmed diagnosis of GD. Detailed study protocol has been described previously.(XX)

3.2. Feature Extraction, Model Training, and Validation

Supplementary Figure S1 shows the convergence of the MSE when training the autoencoders for feature extraction. Recordings from each ECG lead were separately used for training and for each lead the convergence was quick, taking a minimum of 17 (lead aVF) to a maximum of 156 (lead aVR) epochs (Supplementary Table S1). The overall estimate of MSE across the leads ranged from a minimum of 5.74x10^-4 (lead aVL) to a maximum of 1.80x10^-3 (lead V5). When we used the features extracted by the encoders and basic as well as HRV features (total number of features = 818) to train a random forest model on a SMOTE-balanced dataset in a 10-fold cross-validation framework on the training set only, we found that the model performed consistently highly across the folds (Supplementary Figure S2) and the best performing hyperparameter combination was: maximum tree depth = 4, maximum number of features = 3, minimum samples per leaf = 1, minimum samples per split = 2 and number of estimators = 75.

The classification performance (point and 95% CI estimate) of this model averaged across the 10-folds was as follows: AUC, 0.9198 (0.9164 - 0.9233), accuracy, 0.8188 (0.8139 - 0.8237); precision, 0.7370 (0.7316 - 0.7424); recall, 0.9921 (0.9900 - 0.9942); and F1-ratio, 0.8457 (0.8422 - 0.8492). For validation, we used the held-out test set (N = 1,370 with 675 instances of GD) and observed that the model provided comparable and highly accurate prediction of GD status. The classification performance (point and 95% CI estimate) of this model in this held-out test set was as follows: AUC, 0.9152 (0.9005 – 0.9299), accuracy, 0.8445 (0.8244 – 0.8627); sensitivity, 0.8904 (0.8646 – 0.9118); specificity, 0.8000 (0.7656 – 0.8281); precision, 0.8122 (0.7624 – 0.8387); and F1-ratio, 0.8495 (0.8295 – 0.8495) (Figure 3A and 3B). For further analyses, we dubbed this model as the RF-ECG model.

As a final step in validation, we used the model trained on the synthetic data to predict the probability of GD in the original dataset (N = 1,036; 123 instances of GD). In this dataset, the RF-ECG model yielded following metrics (point and 95% confidence interval) of classification performance: AUC, 0.8989 (0.8765 – 0.9213); accuracy 0.7828 (0.7567 – 0.8069); sensitivity 0.8780 (0.8085 – 0.9247); specificity 0.7700 (0.7416 – 0.7961); precision 0.3396 (0.2898 – 0.3933) and F1-ratio 0.4898 (0.4309 – 0.5455). We posited that the decrease in precision and F1-ratio observed in this dataset compared to the held-out test set in the augmented data was primarily due to the stark difference in the prevalence of GD in the two datasets. To that end, we estimated the likelihood ratio of a positive test (LR+) and negative test (LR-) – metrics that are known to be prevalence-agnostic (XX). We observed that in this original dataset, the LR+ was 3.82 (95% CI 3.33 – 4.37) and the LR- was 0.16 (0.10 – 0.25). These estimates were comparable to the LR+ and LR- for the RF-ECG model in the held-out test set of the augmented data: 4.45 (95% CI 3.82 3 – 5.18) and 0.14 (0.11 – 0.17), respectively.

3.3. Comparison of Our Model with Existing Models of GD Risk Prediction

We first investigated if the classification performance of the RF-model improved further by adding statistically significant baseline characteristics (shown in Table 1). The results of these analyses are shown in Figure 3E. The model based on baseline characteristics alone (the MERGD model, blue colored curve in Figure 3E) had a moderate AUC of 0.6670 (95% CI 0.6118 – 0.7221), whereas the AUC (95% CI) for the RF-ECG model alone in the original dataset (red colored curve in Figure 3E) was 0.8989 (0.8765 – 0.9213). However, when this model was combined with the predictions from the MERGD model (green colored curve in Figure 3E), the estimated AUC (95% CI) improved to 0.9103 (0.8889 – 0.9318). This improved AUC (total increase in AUC = 0.0114) was statistically significant (p = 0.0298).

Next, we compared the GD prediction performance of the RF model to that of other externally validated, established risk scoring systems. Specifically, we included the models by Naylor et al [30], the Monash model [7], the van Leeuwen model [29] and the Nanda et al [30] models along with the MERGD and the RF-ECG models alluded to above. The results of these analyses (Figure 3F) showed that all the models based on clinical and socio-demographic characteristics had a comparable AUC (ranging from a minimum of 0.5986 for the Monash model to a maximum of 0.6670 for the MERGD model). However, the classification performance of the RF-ECG model far exceeded that of any of the remaining five models.

3.4. Feature Importances and Model Explanation

In a quest to understand the reasoning behind the predictions provided by the RF-ECG model, we found that the elbow for the plot of feature importances and feature ranks was found at a rank of 40 (4.88%, Figure 4A). Therefore, we focused on the top 40 most important features (which cumulatively accounted for an importance of 28.89%) of all the features. Full list of all the features and their importances (along with their distribution in the dataset) is provided in Supplementary Table S2.

As shown in Figure 4B, an overview of the important features indicated that of the top 40 features, three (mean heart rate, QT interval and P-wave axis) represented basic ECG features reported by the recording device, seven features were related to HRV (six from time domain and one from frequency domain) while the remaining 30 features were identified by the autoencoders. Of these 30 AE-related features (Figure 4C), majority (12, 40%) were in the lead V2 and lead V4 (8, 26.67%).

Figure 4D shows the distribution of the comparative feature values in women who did or did not develop GD. The log percent difference (between GD- and no-GD participants) showed that women with GD had a significantly higher heart rate and a rightward shift of the P-wave axis. The short QT interval observed in women with GD was consistent with a higher heart rate (features corroborating findings shown in Table 1). The six time-domain HRV-related features were PRC80NN (80th percentile of all normal heartbeat-to-heartbeat intervals), MedianNN (median of all normal heartbeat-to-heartbeat intervals), PRC20NN (20^th percentile of all normal heartbeat-to-heartbeat intervals), MaxNN (maximum of all normal heartbeat-to-heartbeat intervals), PNN20 (proportion or of successive heartbeat intervals that differ by more than 20 milliseconds), and TINN (triangular interpolation of the NN interval histogram). The only frequency-domain feature that was included in the top 40 features was VLF (very-low frequency power). As shown in Figure 4D, all the HRV-related features showed lower values in the GD patients as compared to the non-GD participants, indicating reduced HRV associated with GD. Together, the feature importance analysis indicated that patients with GD had higher heart rate, rightward deviation of the P-wave axis, reduced heart rate variability and a concentration of convolutionally detected features on the V2 and V4 leads.

4. Discussion

We developed and validated a deep learning/machine learning model to predict incident GD early during pregnancy using electrocardiography. The model predicted GD with high degree of accuracy and outperformed the existing clinical risk stratification systems for GD. While electrocardiographic changes have been consistently shown to be associated with the risk of GD, to our knowledge, this is first such model that taps subtle ECG features to predict the risk of GD 10.41 weeks ahead of a confirmatory OGTT. This lead time can be efficiently used in clinical practice to meliorate the risk of incident GD. Lastly, our results also point towards the possibility of combining the ECG-based model with other clinical and socio-demographic variables to further improve its classification performance.

Zhang et al[33] conducted an elegant meta-analysis of 25 studies using machine learning models for prediction of gestational diabetes. These models used tabular, clinical data extracted from electronic health records containing combinations of the following variables: maternal age, family history of diabetes, fasting plasma glucose, pre-pregnancy BMI, history of diabetes, serum triglycerides, hemoglobin A1c, systolic blood pressure and high-sensitivity C-reactive protein. The pooled AUC for these machine learning models was reported to be 0.84 but the dichotomized performance for GDM diagnosis had a sensitivity of 0.69 and specificity of 0.75. However, currently there are no studies that use digital electrocardiographic input and combine these with artificial intelligence for prediction of GD. Moreover, the sensitivity and specificity of our model was superior to the estimates reported by Zhang et al.[33]

The biological plausibility of ECG as a marker of GD rides on previous studies. Pregnancy offers a high-stress environment to the mother’s heart.[34] As a result, the maternal heart is in a process of adaptation throughout pregnancy. Thirunavukarasu et al[35] demonstrated that a higher left ventricular (LV) mass, lower LV end-diastolic volume and lower global longitudinal shortening are characteristic of patients with GD. These changes are suggestive concentric LV hypertrophy consistent with hypertensive heart disease. Heart rate variability measures (QT dispersion)[14,36] and ECG-based body surface mapping[37] have also been shown to be associated with GD. Gasic et al [38] have also shown that features of cardiac autonomic neuropathy (CAN) are more frequent in GD than in non-GD women. Notably, GD is often associated with insulin resistance which exercises a metabolic strain on the heart muscle in addition to the preload due to plasma volume expansion.[34] In this context, two findings from our study deserve a mention. First, we found features suggestive of an increased right heart preload in GD (P-wave axis shift and increased heart rate). This finding is concomitant with the expected expansion of plasma early during pregnancy.[39] Second, we observed a reduced HRV with reduced very-low frequency power in women with GD – a finding that corroborates a likely autonomic axis to the involvement of maternal heart in GD.[40]

There are some important limitations of the study that preclude immediate generalization of the findings. First, despite the robustness of cross validation and of using the original dataset as an external dataset (Figure 1), the real test of the algorithm will be in external, independent datasets. Unfortunately, such well-curated datasets that include digitally collected ECGS are currently unavailable. Till the time the models are externally validated in independent data sets, generalization of the model cannot be done. Second, the autoencoder-related features are essentially uncharacterized. Therefore, a full understanding of the ECG features contributing to GD is still not possible at this time. Further studies are needed (e.g. the morphing methods used by Obermeyer et al [41]) to innovatively characterize the subtle ECG features that are associated with GD. Third, GD is mixed bag od varying definitions and thresholds. The ground truth definition used in this study was based on a previous paper focusing on the study population. This definition of GD is not universal and therefore the model will need to be made adaptable to the varying definitions of GD. At this time, our study is a proof-of-principle endeavor to demonstrate the utility of deep learning / machine learning in developing an ECG-based biomarker for GD.

5. Conclusions

Akin to the fact that a GCT is a gate-keeper first step recommended [2] in the diagnosis of GD and that the GCT can also be informative in the context of future cardiovascular health,[42] we have developed a non-invasive, early and accurate method to predict the likelihood of GD that promises substantial lead time. The model detected features consistent with clinical expectations and more accurately than existing risk stratification systems for GD. If the classification performance of our model is supported by external datasets in future studies, then our approach has a potential to provide a feasible alternative strategy for GD diagnosis, especially in resource limited settings.