The level of glycated hemoglobin (HbA_1c) is used to measure the average glucose concentration in red blood cells [1,2]. Unlike other glucose blood tests, such as random blood sugar (RBS) and fasting blood sugar (FBS), HbA_1c provides a long-term measure of a patient’s blood glucose levels [3]. The HbA_1c test can therefore provide physicians with a reliable means of monitoring a patient’s hyperglycemia without requiring the patient to undertake overnight fasting prior to being tested.

A concentration of 6.5% for the HbA_1c in patient blood is considered as the cutoff point for the diagnosis of diabetes [4]. However, patients with a concentration of less than 6.5% are not completely excluded from a diabetes diagnosis, as the range of elevation levels (5.7%≤ HbA_1c <6.5%) can indicate the future onset of diabetes. Therefore, HbA_1c can act as an early predictor for the potential development of type-2 diabetes mellitus (T2DM) [2]. Ackermann et al [3] suggested using the HbA_1c test as a measure for identifying those adults who are at a greater risk of developing T2DM in the future.

Research has shown that reducing HbA_1c levels can significantly reduce the possibility of developing serious complications. Hence, close monitoring of HbA_1c levels is recommended for all diabetic patients and those with the potential for developing diabetes [5]. It is also suggested that diabetic and nondiabetic patients with raised HbA_1c levels should be clinically checked and monitored as a preventive intervention to avoid developing T2DM [6].

Currently, the clinical data collected from patient visits consists of a set of readings for vital signs and lab tests, diagnoses, physicians’ notes, and treatments that are stored in electronic health records (EHRs). These are collected on an irregular basis, according to clinical needs, and stored with an associated time stamp.

In recent years, machine learning models have shown powerful capabilities for analyzing and understanding complex data across a wide variety of applications. Our research question for this study was as follows: “Can HbA_1c prediction be improved by using machine learning with longitudinal data that are normally available in EHR systems?”

This paper reports an investigation into the performance of machine learning models to predict current HbA_1c levels as a binary classification problem using EHR data. Nondiabetic patients with an HbA_1c level of 5.7% or more are considered to have an elevated HbA_1c, while those with levels lower than this are considered normal. The models combine current visit data with extra features (independent variables) extracted from previous visits by patients. We used explainable methods to rank the features in order of their importance to the decision made by each of the models. To the best of our knowledge, this study is the first to employ machine learning models that use longitudinal data from EHR systems for the purpose of HbA_1c elevation risk prediction. This study is also the first to use explainable machine learning techniques to explain the classification decisions made by black box models, support vector machine (SVM), and multilayer perceptron (MLP), in predicting HbA_1c elevation risk (≥5.7%), in order to better understand the behavior of the model.

EHR data have been intensively investigated for a variety of medical decision support tasks [7]. These tasks include the analysis of complex patterns and prediction of major medical events (for example, diagnostic imaging and gene interactions) [8,9]. Several studies have demonstrated the successful employment of EHR data with prediction models [10]. For instance, machine learning has been intensively used with EHR data in diagnosing diabetes and discovering its related patterns [11-15]. However, we are not aware of any studies that have explored machine learning models for the prediction of current elevated HbA_1c levels using EHR data from a nondiabetic population or the impact of patient longitudinal data on the effectiveness of such predictive machine learning models.

Several studies have investigated the association between HbA_1c levels and clinical variables using statistical models [16,17]. A study by Rose et al [18] discussed the correlation between RBS and HbA_1c levels. Stanley et al [19] used a linear regression model for imputation of missing HbA_1c data. Their model calculates HbA_1c levels for patient records with missing HbA_1c values as continuous and categorical values and uses 4 predictors extracted from an EHR system—RBS, FBS, age, and gender—as predictors to calculate the level of HbA_1c for a diabetic population. Simone et al [20] used linear regression models to predict HbA_1c levels after 6 years for nondiabetic patients using different populations.

A study by Wells et al [21] in 2018 was the first to focus on predicting current HbA_1c elevation levels for nondiabetic patients through use of an EHR data set. Multiple logistic regression (MLR) was employed to calculate the probability of a patient having an elevated HbA_1c level (≥5.7%). The data set was extracted from an EHR system used in the United States. The authors used 8 independent variables fitted to the model using restricted cubic splines with 3 knots to formulate the final equation. The performance of the MLR model was compared to that of the models used by Baan et al [22] and Griffin et al [23]. However, the models by Baan and Griffin aimed at predicting the onset of patients’ diabetes rather than predicting HbA_1c levels for nondiabetic patients. In addition, the experimental data set used by Wells et al to train and test their model was imbalanced with 74% of the samples having normal HbA_1c levels (5.7%) and only 26% of the samples having elevated HbA_1c levels (≥5.7%).

We performed a differentiated replication of the study by Wells et al [21] using the more balanced King Abdullah International Medical Research Center (KAIMRC) data set [24]. Although the significant variables identified in our replication were in general agreement with those of the original study, there were some differences in the ranking of importance for these, suggesting that such models do need to be “tuned” to the characteristics of different populations.

To study the impact of using advanced predictive models with EHR data to predict current HbA_1c levels, we employed the MLR, random forest (RF), SVM, and logistic regression (LR) models, as well as a deep learning model, MLP [25]. The problem was formulated into a binary classification problem whereby the target variable, HbA_1c level, was encoded as 1 when the level of HbA_1c was 5.7% or more and with 0 otherwise. The results obtained from using these models were compared to those obtained from employing the model used by Wells et al with the KAIMRC data set (detailed in the Data Set subsection).

The performance of the models was investigated using current visit data only and with additional longitudinal data from current and previous visits. The performance of each model was evaluated using measures commonly employed in clinical applications. For the SVM and MLP models, the relative importance of the features was also calculated using explainable machine learning techniques.

Using black box machine learning models in health care can have adverse effects on the trust and confidence placed in their outcomes; the risk of misclassification is potentially too high for clinicians to confidently use black box models for high risk health care decisions, and not being able to interpret a model’s decision exacerbates this problem [26]. Explainable methods for machine learning models allow interpretable outcomes that can expose the reasons behind the decision made by the model [27]. This transparency provides both health professionals and patients with the confidence and trust in the outcome of the models. The widely used Shapley Additive Explanations (SHAP) values [28] and local interpretable model-agnostic explanations (LIME) score [29] techniques have therefore been used to provide a degree of transparency to our deep learning model.

SHAP values are derived from Shapley values used in game theory and provide a method of calculating the contribution of each feature (variable) to the final prediction via the GradientSHAP approximation. This is achieved for each feature by comparing the prediction the model makes when the feature is present with the prediction obtained when the feature takes some baseline value [28]. Consequently, the SHAP values for a given input “explain” how each feature affects the output of the model when compared to the baseline (or “default”) output of the model. We used SHAP values to interpret our black box models, so they could be efficiently calculated, and their use enabled a global view of the model to be constructed through the computation of SHAP values from across the whole data set.

SHAP values were computed using the feature’s mean marginal contribution across different coalitions of all features. SHAP values themselves are computationally intensive to compute, and so approximation methods are commonly used when calculating the values.

To ensure that the SHAP values we calculated were not too greatly affected by the approximation method used, we also computed the LIME [29] scores for the models across the entire data set. LIME tries to estimate locally faithful linear explanations (ie, explanations that correspond to how the model behaves around the instance being explained) for any classifier. LIME achieves this by creating local linear classifiers that approximate the behavior of the original model in the vicinity of the data being explained. As linear models are inherently interpretable through their parameters, they can be used to generate explanations of the original model. Both SHAP and LIME have the advantage that they are model-agnostic techniques, and so we were able to apply both methods to both of our black box classification models (SVM and MLP).

The data used in this study were taken from the KAIMRC data set. The data were collected from King Abdulaziz Medical City located in the central and western regions of Saudi Arabia, an area which has been ranked second in the Middle East and seventeenth in world in diabetes prevalence by the World Health Organization (WHO) [30]. According to the International Diabetes Federation, the diabetes prevalence rate in Saudi Arabia is 18.3%. Therefore, the availability of the data from this population provides considerable opportunities for research into the early prediction of diabetes.

The data set contains a full history of patient details, vital signs, and lab test readings for each patient visit for the period from 2016 to the end of 2018. As the aim of this study was to identify nondiabetic patients that are at a high risk of HbA_1c elevation, all patients previously diagnosed with hyperglycemia were excluded from the experimental data set. The remaining cohort formed our experimental data set and was categorized by using the American Diabetes Association’s guidelines [31], in which patients with HbA_1c readings of more than 5.7% are considered as being in the prediabetic range, while those with less than 5.7% are considered to be in the normal range.

Most medical data sets are imbalanced [32-34]. These imbalances occur when the proportion of one class of patients in the data set is greater than its counterpart class [35,36]. However, unusually, our experimental data set was not imbalanced. Slightly over half of the patients in our experimental data set (9826/18,844, 52.14%) were found to have elevated levels of HbA_1c (≥5.7%) while 47.86% (9018/18,844) of patients had normal HbA_1c levels (<5.7%). This can be ascribed to the high incidence of diabetes in the region from which the data set was collected [37].

A detailed illustration of the patients’ class distribution (HbA_1c levels) by age groups and gender is shown in Figure 1. This shows that as the age of patients increased, so did the proportion of patients who had elevated HbA_1c levels. The data set also exhibited a balanced gender distribution, with 49.40% (9308/18,844) of the patients being male and 50.60% (9536/18,844) being female. However, the proportion of male patients with elevated levels of HbA_1c (≥5.7%) was greater than that of the female patients. Also, female patients with normal levels of HbA_1c (<5.7%) made more visits than did males. Table 1 shows the profile for the distribution of HbA_1c elevation levels organized by gender.

Six main variables (features) were extracted from the KAIMRC EHR data set to be used in this study. These features, which were selected first for their theoretical association with hyperglycemia and second for their availability in the KAIMRC data set, were the following: age, BMI, estimated glomerular filtration rate (eGFR), RBS, total cholesterol, and non–high-density lipoprotein. The lab codes of the features used are available in Multimedia Appendix 1 Table S1. The descriptive statistics (using the data for the current visit only for unique patients), units, and P values for the selected features are presented in Table 2.

It is very common in clinical practice that physicians may require some lab tests and vital signs to be frequently recorded. In these cases, the average value of all readings taken on a given day (the basic time interval used for this study) was used. For inpatient visits, only data for the first day were considered, and, where there were missing values, the first available values from the visit were used.

For the purpose of this study, we aimed at predicting the HbA_1c levels (≥5.7%) for current (last) patient visits only. Unlike the sampling approach used by Wells et al, which was based on independent hospital visits for patients (including for the same patients), the sampling approach used in this study included independent patients to ensure only unseen patients data were used for testing the models. Although we aimed to identify patients with elevated levels of HbA_1c from a nondiabetic population, patients previously diagnosed with diabetes were excluded. We also excluded nonadult patients and those with erroneous or missing values [24]. Figure 2 shows the details of the tasks performed to refine the sample selection. This resulted in a reduction in the size of the experimental data set from 114,057 patients with 750,709 visits to 18,844 unique patients with 157,600 visits.

The inputs (input features space) for the models used in this study were continuous values. Values for age, eGFR, RBS and total cholesterol features were directly available in the KAIMRC data set. The values for the BMI and non–high-density lipoprotein variables were calculated from other available features using the formulae in Multimedia Appendix 2.

The input structure for the deep learning model was organized as a matrix, based on current and previous time-stamped patient visits. It contained the current visit data concatenated with approximated values for the selected features from all previous visits, which we refer to as the “Approximated Time Series Data”.

Each patient visit was described by the selected features, represented as x₁, x₂ …, x_n. These features were formed as episodes based on the time-stamped values available in each visit (v_i).

As a baseline comparison, we employed the MLR model used by Wells et al [21], and compared the results from this with those from 4 commonly used machine learning models.

The MLR model is used to create a mathematical equation that can best calculate the probability of a value by assigning weights (coefficients) to the independent variables (features) based on their importance [40]. In this study we employed the same approach used by Wells et al by which the continuous features were fitted into the MLR model using restricted cubic splines technique with 3 knots. When we used the longitudinal input, the variables that caused collinearity were excluded.

Random forest is an algorithm very commonly used for classification. It combines several decision trees that are generated during the training process. Each decision tree is trained using a random subset of the training data set. The final classification is then based on the majority voting results of all generated decision trees [41]. The quality function used in the employed RF model is the Gini importance, with a value of 100 for the number of tree parameters.

Logistic regression is commonly used to solve binary classification problems. It calculates the odds ratio of the variables and is similar to MLR but uses a binomial distribution of the dependent variable (ie, more than 1). Thus, it includes a logit function that handles different types of relationships between the dependent and independent variables [42,43].

Support vector machine was introduced by Vapnik [44] in 1998. It can solve both classification and regression problems. It uses the training feature space to decide on the separation boundaries (hyperplane) that best divides the training data set into regions, 1 for each class. The very close points to the hyperplanes are the support vectors. SVMs also use kernels to help enhance class separation by mapping the training features into a higher dimensional space with an increased number of dimensions [44,45]. The kernel function used in the SVM model employed is a radial base function with a value of 1 for the cost parameter (C).

A multilayer perceptron, also known as a feed-forward neural network, is one of the most common deep learning approaches. It is mainly used to address supervised learning problems by learning the dependencies between the input layer (the features or variables) and output layer (the classification decision) using a fully connected hidden layer in between. The layers, including hidden ones, contain a number of neurons that are connected to the neurons of the next and previous layers via weights and nonlinear functions. MLP uses a backpropagation algorithm to update the weights and biases within the hidden layers to minimize the output error rate [25,46].

To optimize the MLP model, fine-tuning of the structure and hyperparameters was performed and involved the number of hidden layers and neurons, activation functions, optimizers, and loss functions. The optimized structure of the MLP model used in this study contained 3 hidden layers. The number of neurons in the hidden layers were 48, 48, and 24, respectively. The final layer (the output layer) contained 2 neurons for the final output of the model (Y1 for normal HbA_1c or Y2 for elevated HbA_1c). A rectified linear unit activation function was used in the 3 hidden layers, while a sigmoid was used in the output layer. The detailed structure of the MLP model is shown in Figure 3. The model was trained using an Adam optimizer with mean squared error as the loss function.

The models all employed the same data preprocessing, training, and testing techniques. The models were validated using the 10-fold cross-validation technique. The k-fold cross-validation is one of the most commonly used approximation approaches for validating the obtained results [47,48]. For the MLP model, 100 epochs were used to train each fold.

As our measure for evaluating and comparing the performance of the proposed models, we used the area under the receiver operating characteristic (AUC-ROC) curve, which is equal to the concordance statistic [49]. We also report values for a set of measures that are commonly used in clinical applications: balanced accuracy (that calculates the recall average for each class), overall accuracy, F score, precision, and precision-recall area under the curve (PR-AUC).

To determine the importance that the black box models (SVM and MLP) place upon each variable, we first computed the SHAP values and LIME scores for all samples in our data set and then calculated the average absolute SHAP value and LIME score for each predictor.

EHR systems were adopted for the purpose of improving health care outcomes and were not originally intended for research purposes [19]. Patient data stored in EHR systems can be obtained at irregular intervals, as lab instructions are carried out with different frequencies based on the physician's decisions and a patient’s visit patterns. It is very common that medical data extracted from EHR systems suffer from problems such as irregularity, incompleteness, and noisy and imbalanced data [13]. These can be challenging obstacles for any technology used for predictive analytics.

In our study, the sampling approach used did not affect the balanced nature of the data set used. As shown in Figure 2, there were 56,185 unique patients present before removal of the records with 1 or more missing values. The number of unique patients with elevated HbA_1c levels (≥5.7%) before removal of the incomplete records was 27,354, resulting in a retention of 48.68% (27,354/56,185). The number of unique patients with normal HbA_1c levels was 28,831, resulting in a retention of 51.32% (28,831/56,185). We would argue that the absence or the presence of the HbA_1c readings is not random, as the sample was collected from the population of Saudi Arabia and thus the likelihood of a patient taking an HbA_1c test is large because of the prevalence of diabetes in this country [51]. This may affect the reproducibility of this work using different populations from different countries especially those with lower rates of diabetes.

It is hoped that these outcomes will encourage further investigation into the predictability of current HbA_1c levels (≥5.7%) using more of the readings normally provided in EHR data. For example, other important readings such as FBS and triglycerides have shown clinical correlations with diabetes [52]. In addition, our data set contained only 3 years of patient data, which limits the number of patient visits recorded. Figure 11 shows the number of visits made by patients from 2016 to 2018, while Figure 12 details the number of visits made by patients (after removal of the outliers) over HbA_1c levels. Both figures show that the majority of the patients have made relatively few visits: 52% (8713/16818) of the patients made 4 visits or fewer over the 3 years (1.3 visit per year). This also justifies the size of the sliding window (s=3) as the optimal input size for the models used. However, we hypothesize that the longitudinal behavior of the features used can be enriched by including more values obtained over longer periods. Therefore, incorporating more features and their longitudinal behavior over longer periods into the models used in this study would likely improve the prediction performance of our chosen models.

Variations in the data or model produce slightly different attribution values. However, due to the critical nature of many health care applications, it is always important to verify that the models make “sensible” predictions. Without the use of SHAP/LIME, this would be hard to verify for any nonlinear model. Although it is possible to see that the models have high performance, we would be unable to verify that a model is not making spurious correlations. Furthermore, through the use of SHAP, we can verify that MLPs trained on the longitudinal data are learning to use the extra information contained in the longitudinal data (as indicated by the higher importance of eGFR), allowing us to pinpoint the reason these models gain higher performance.

To investigate the effect of temporal dependencies in the data, this study investigated the use of other deep learning models along with the MLP, including long short-term memory (LSTM) and bidirectional LSTM [25,53] for HbA_1c prediction. Table 5 reports the results of using these models. The MLP model achieved similar performance to the LSTM and bidirectional LSTM models according to all reported measures. This suggests that directly modeling the temporal dynamics in the data is not very helpful. This could be due to the short lengths of the time series or a too-weak temporal dependency.

Generalizing our findings using other data sets is challenging because of the accessibility and privacy restrictions that apply to medical data sets. For this reason, and because of the lack of similar studies that have used machine learning for HbA_1c prediction with EHR data, comparing the performance achieved by the models outlined in this study with those developed by other researchers will require the availability of alternative anonymized data sets.

We believe that this study is the first to investigate the performance of machine learning models used with EHR data for predicting current HbA_1c elevation risk (≥5.7%) for nondiabetic patients. It is also the first to investigate employing the longitudinal data that are normally stored on EHR systems to enhance the prediction of HbA_1c elevation levels. Our findings show that the MLP model achieves better results when a patient’s longitudinal data are combined with current visit data, and the use of longitudinal data also affects the relative importance for the predictors used.

As this work formed a continuation of previous work [24], we avoided changing the sampling approach used. However, studying the impact of applying different sampling approaches could be valuable to explore in future work as would the use of a larger data set with more variables and the recording of longitudinal behavior over longer periods.