A function approximation approach to the prediction of blood glucose levels

The problem of real time prediction of blood glucose (BG) levels based on the readings from a continuous glucose monitoring (CGM) device is a problem of great importance in diabetes care, and therefore, has attracted a lot of research in recent years, especially based on machine learning. An accurate prediction with a 30, 60, or 90 minute prediction horizon has the potential of saving millions of dollars in emergency care costs. In this paper, we treat the problem as one of function approximation, where the value of the BG level at time $t+h$ (where $h$ the prediction horizon) is considered to be an unknown function of $d$ readings prior to the time $t$. This unknown function may be supported in particular on some unknown submanifold of the $d$-dimensional Euclidean space. While manifold learning is classically done in a semi-supervised setting, where the entire data has to be known in advance, we use recent ideas to achieve an accurate function approximation in a supervised setting; i.e., construct a model for the target function. We use the state-of-the-art clinically relevant PRED-EGA grid to evaluate our results, and demonstrate that for a real life dataset, our method performs better than a standard deep network, especially in hypoglycemic and hyperglycemic regimes. One noteworthy aspect of this work is that the training data and test data may come from different distributions.


Introduction
Diabetes is a major disease affecting humanity. According to the 2020 National Diabetes Statistics report [3], more than 34,000,000 people had diabetes in the United States alone in 2018, contributing to nearly 270,000 deaths and costing nearly 327 billion dollars in 2017. There is so far no cure for diabetes, and one of the important ingredients in keeping it in control is to monitor blood glucose levels regularly. Continuous glucose monitoring (CGM) devices are used increasingly for this purpose. These devices typically record blood sugar readings at 5 minute intervals, resulting in a time series. This paper aims at predicting the level and rate of change of the sugar levels at a medium range prediction horizon; i.e., looking ahead for 30-60 minutes. A clinically reliable prediction of this nature is extremely useful in conjunction with other communication devices in avoiding unnecessary costs and dangers. For example, hypoglycemia may lead to loss of consciousness, confusion or seizures [9]. However, since the onset of hypoglycemic periods are often silent (without indicating symptoms), it can be very hard for a human to predict in advance. If these predictions are getting communicated to a hospital by some wearable communication device linked with the CGM device, then a nurse could alert the patient if the sugar levels are expected to reach dangerous levels, so that the patient could drive to the hospital if necessary, avoiding a dangerous situation, not to mention an expensive call for an ambulance.
Naturally, there are many papers dealing with prediction of blood glucose (BG) based not just on CGM recordings but also other factors such as meals, exercise, etc. In particular, machine learning has been used extensively for this purpose. We refer to [9] for a recent review of BG prediction methods based on machine learning, with specific focus on predicting and detecting hypoglycemia.
The fundamental problem of machine learning is the following. We have (training) data of the form {(x j , f (x j )+ j )}, where, for some integer d ≥ 1, x j ∈ R d are realizations of a random variable, f is an unknown real valued function, and the j 's are realizations of a mean zero random variable. The distributions of both the random variables are not known. The problem is to approximate f with a suitable model, especially for the points which are not in the training data. For example, in the prediction of a time series x 0 , x 1 , . . . with a prediction horizon of h based on d past readings can be (and is in this paper) viewed as a problem of approximating an unknown function f such that x j+h = f (x j−d+1 , . . . , x j ) for each j. Thus, x j = (x j−d+1 , . . . , x j ) and y j = x j+h . Such problems have been studied in approximation theory for more than a century. Unfortunately, classical approximation theory techniques do not apply directly in this context, mainly because the data does not necessarily become dense on any known domain such as a cube, sphere, etc., and we cannot control the locations of the points x j . This is already clear from the example given above -it is impossible to require that the BG readings be at a pre-chosen set of points.
Manifold learning tries to ameliorate this problem by assuming that the data lies on some unknown manifold, and developing techniques to learn various quantities related to the manifold, in particular, the eigen-decomposition of the Laplace-Beltrami operator. Since these objects must be learnt from the data, these techniques are applicable in the semi-supervised setting, in which we have all the data points x j available in advance. In [8], we have used deep manifold learning to predict BG levels using CGM data. Being in the semi-supervised setting, this approach cannot be used directly for prediction based on data that was not included in the original data set we worked with.
Recently, we have discovered [7] a more direct method to approximate functions on unknown manifolds without trying to learn anything about the manifold itself, and giving an approximation with theoretically guaranteed errors on the entire manifold; not just the points in the original data set. The purpose of this paper is to use this new tool for BG prediction based on CGM readings which can be used on patients not in the same data sets, and for readings for the patients in the data set which are not included among the training data.
In this context, the numerical accuracy of the prediction alone is not clinically relevant. The Clarke Error-Grid Analysis (C-EGA) [2] is devised to predict whether the prediction is reasonably accurate, or erroneous without serious consequences (clinically uncritical), or results in unnecessary treatment (overcorrection), or dangerous (either because it fails to identify a rise or drop in BG or because it confuses a rise in BG for a drop or vice versa). The Prediction Error-Grid Analysis (PRED-EGA), on the other hand, categorizes a prediction as accurate, erroneous with benign (not serious) consequences or erroneous with potentially dangerous consequences in the hypoglycemic (low BG), euglycemic (normal BG), and hyperglycemic (high BG) ranges separately. This stratification is of great importance because consequences caused by a prediction error in the hypoglycemic range are very different from ones in the euglycemic or the hyperglycemic range. In addition, the categorization is based not only on reference BG estimates paired with the BG estimates predicted for the same moments (as C-EGA does), but also on reference BG directions and rates of change paired with the corresponding estimates predicted for the same moments. It is argued in [14] that this is a much better clinically meaningful assessment tool than C-EGA, and is therefore the assessment tool we choose to use in this paper.
Since different papers on the subject use different criterion for assessment it is impossible for us to compare our results with those in most of these papers. Besides, a major objective of this paper is to illustrate the use and utility of the theory in [7] from the point of view of machine learning. So, we will demonstrate the superiority of our results in the hypoglycemic and hyperglycemic regimes over those obtained by a blind training of a deep network using the MATLAB Deep Learning Toolbox. We will also compare our results with those obtained by the techniques used in [10,6,8]. However, we note that these methods are not directly comparable to ours. In [10], the authors use the data for only one patient as training data, which leads to a meta-learning model that can be used serendipitously on other patients, although a little bit of further training is still required for each patient to apply this model. In contrast, we require a more extensive "training data", but there is no training in the classical sense, and the application to new patients consists of a simple matrix vector multiplication. The method in [6] is a linear extrapolation method, where the coefficients are learnt separately on each patient. It does not result in a model that can be trained once for all and applied to other patients in a simple manner. As pointed out earlier, the method in [8] requires that the entire data set be available even before the method can start to work. So, unlike our method, it cannot be trained on a part of the data set and applied to a completely different data set, not available at the beginning.
We will explain the problem and the evaluation criterion in Section 3. Prior work on this problem is reviewed briefly in Section 2. The methodology and algorithm used in this paper are described in Section 4. The results are discussed in Section 5. The mathematical background behind the methods described in Section 4 is summarized in Appendix A.

Prior work
Given the importance of the problem, many researchers have worked on it in several directions. Below we highlight the main trends relevant to our work.
The majority of machine learning models (30% of those studied in the survey paper [9]) are based on artificial neural networks such as convolutional neural networks, recurrent neural networks and deep learning techniques. For example, in [17], the authors develop a deep learning model based on a dilated recurrent neural network to achieve 30-minute blood glucose prediction. The authors assess the accuracy of their method by calculating the root mean squared error (RMSE) and mean absolute relative difference (MARD). While they achieve good results, a limitation of the method is the reliance on various input parameters such as meal intake and insulin dose which are often recorded manually and might therefore be inaccurate. In [11], on the other hand, a feed-forward neural network is designed with eleven neurons in the input layer (corresponding to variables such as CGM data, the rate of change of glucose levels, meal intake and insulin dosage), and nine neurons with hyperbolic tangent transfer function in the hidden layer. The network was trained with the use of data from 17 patients and tested on data from 10 other patients for a 75-minute prediction, and evaluated using C-EGA. Although good results are achieved in the C-EGA grid in this paper, a limitation of the method is again the large amount of additional input information necessary to design the model, as described above.
A second BG prediction class employs decision trees. The authors of [13], for example, use random forests to perform a 30-minute prediction to specifically detect postprandial hypoglycemia, that is, hypoglycemia occurring after a meal. They use statistical analyses like sensitivity and specificity to evaluate their method. In [16], prediction is again achieved by random forests and incorporating a wide variety of additional input features such as physiological and physical activity parameters. Mean Absolute Percentage Error (MAPE) is used as assessment methodology. The authors of [5] use a classification and regression tree (CART) to perform a 15-minute prediction using BG data only, again using sensitivity and specificity to test the performance of their method.
A third class of methods employ kernel-based regularization techniques to achieve prediction (for example, [10] and references therein), where Tikhonov regularization is used to find the best least square fit to the data, assuming the minimizer belongs to a reproducing kernel Hilbert space (RKHS). Of course, these methods are quite sensitive to the choice of kernel and regularization parameters. Therefore, the authors in [10] develop methods to choose both the kernel and regularization parameter adaptively, or through meta-learning ("learning to learn") approaches. C-EGA and PRED-EGA are used as performance metrics.
Time series techniques are also employed in the literature. In [12], a tenth-order auto-regression (AR) model is developed, where the AR coefficients are determined through a regularized least square method. The model is trained patient-by-patient, typically using the first 30% of the patient's BG measurements, for a 30-minute or 60-minute prediction. The method is tested on a time series containing glucose values measured every minute, and evaluation is again done through the C-EGA grid. The authors in [15] develop a first-order AR model, patient-bypatient, with time-varying AR coefficients determined through weighted least squares. Their method is tested on a time series containing glucose values measured every three minutes, and quantified using statistical metrics such as RMSE. As noted in [10], these methods seem to be sensitive to gaps in the input data.
A summary of the methods described above is given in Table 1.
It should be noted that the majority of BG prediction models (63.64% of those analyzed in [9]), including many of the references described above, are dependent on additional input features (such as meal intake, insulin dosage   and physical activity) in addition to historical BG data. Since such data could be inaccurate and hard to obtain, we have intentionally made the design decision in our work to base our method on historical BG data only.

The set-up
We use two different clinical data sets provided by the DirectNet Central Laboratory [1], which lists BG levels of different patients taken at 5-minute intervals with the CGM device; i.e., for each patient p in the patient set P , we are given a time series {s p (t j )}, where s p (t j ) denotes the BG level at time t j . The relevant details of the two data sets are described in Table 2.  Our goal is to predict for each j, the level s p (t j+h ), given readings s p (t j ), . . . , s p (t j−d+1 ) for appropriate values of h and d. We took d = 7 (a sampling horizon t j − t j−d+1 of 30 minutes has been suggested as the optimal one for BG prediction in [6,4]). We tested prediction horizons of 30 minutes (h = 6) (the most common prediction horizon according to the analysis in [9]), 60 minutes (h = 12) and 90 minutes (h = 18).

Data reformulation
For each patient p in the set P of patients in the data sets, the data is in the form of a time series {s p (t j )} Np j=1 of BG levels at time t j , where t j − t j−1 = 5 minutes. We re-organize the data in the form For any patient p ∈ P , we will abbreviate The problem of BG prediction is then seen as the classical problem of machine learning; i.e., to find a functional relationship f such that f (x p,j ) ≈ y p,j , for all p and j = d, · · · , N p − h.

Training data selection
To obtain training data, we form the training set C by randomly selecting (according to a uniform probability distribution) c% of the patients in P . The training data are now defined to be all the data of each patient in C, that is, Since we define the training and test data in terms of patients, choosing the entire data for a patient in the training (respectively, test) data, we may now simplify the notation by writing (x j , y j ) rather than (x p,j , y p,j ), with the understanding that if x j represents a data point for a patient p, y j is the corresponding value for the same patient p. For future reference, we also define

BG range classification
To get accurate results in each BG range, we want to be able to adapt the training data and parameters used for predictions in each range -as noted previously, consequences of prediction errors in the separate BG ranges are very different. To this end, we divide the measurements in C into three clusters C o , C e and C r , based on y j , the BG value at time t j+h , for each x j : In addition, we also classify each x j ∈ Q as follows: Ideally, the classification of each x j ∈ Q should also be done based on the value of y j (as was done for the classification of the training data C). However, since the values y j are only available for the training data and not for the test data, we use the next best available measurement for each x j , namely x j,d , the BG value at time t j .

Prediction
Before making the BG predictions, it is necessary to scale the input readings so that the components of each x j ∈ P are in [−1/2, 1/2]. This is done through the transformation where M := max{x j,k : x j ∈ C} and m := min{x j,k : x j ∈ C}.
WithF n,α =F n,α (Y, x j ) defined as in (A.3) used with training data Y and evaluated at a point x j , we are finally ready to compute the BG prediction for each x j ∈ Q by The construction of the estimatorF n,α involves several technical details regarding the classical Hermite polynomials. For ease of reading, the details of this construction are deferred to the appendix.

Evaluation
To evaluate the performance of the final output f (x j ), x j ∈ Q, we use the PRED-EGA mentioned in Section 3. Specifically, a PRED-EGA grid is constructed by using comparisons of f (x j ) with the reference value y j . This involves comparing f (x j ) with y j as well as the rates of change for all x j ∈ Q. Based on these comparisons, PRED-EGA classifies f (x j ) as being Accurate, Benign or Erroneous.
We repeat the entire process described in Subsections 4.2 -4.5 for a fixed number of trials, after which we report the average of the PRED-EGA grid placements, over all x j ∈ Q and over all trials, as the final evaluation.
A summary of the method is given in Algorithm 1.

Results and discussion
As mentioned in Section 3, we apply our method to data provided by the DirecNet Central Laboratory. Time series for the 25 patients in data set D and the 25 patients in data set J that contain the largest number of BG measurements are considered. These specific data sets were obtained to study the performance of CGM devices in children with Type I diabetes; as such, all of the patients are less than 18 years old. Summary statistics of the two data sets are provided in Table 2. Our method is a general purpose algorithm, where these details do not play any significant role, except in affecting the outcome of the experiments.
We provide results obtained by implementing our method in MATLAB, as described in Algorithm 1. For our implementation, we employ a sampling horizon t j − t j−d+1 of 30 minutes (d = 7), 50% training data (c = 50) (which is comparable to approaches followed in for example [12]) and a total of 100 trials (T = 100). For all the experiments, the function approximation parameters n o , n e and n r referenced in Algorithm 1 were chosen in the range {3, . . . , 7} with α o = α e = α r = 1. We provide results for prediction horizons t j+m − t j of 30 minutes (h = 6), 60 minutes (h = 12) and 90 minutes (h = 18) (again, comparable to prediction horizons in for example [12,13,17]). After testing our method on all 25 patients in each data set separately, the average PRED-EGA scores (in percent) are displayed in Tables 3 and 4. For comparison, Tables 3 and 4 also display predictions on the same data sets using four other methods: (i) MATLAB's built-in Deep Learning Toolbox. We used a network architecture consisting of three fully connected convolutional and output layers and appropriate input, normalization, activation, and averaging layers. As before, we used 50% training data and averaged the experiment over 100 trials, for prediction windows of 30 minutes, 60 minutes and 90 minutes and data sets D and J.
(ii) The diffusion geometry approach, based on the eigen-decomposition of the Laplace-Beltrami operator, we followed in our 2017 paper [8]. As mentioned in Section 1, this approach can be classified as semi-supervised learning, where we need to have all the data points x j available in advance. Because of this requirement, Algorithm 1 Deep Network for BG prediction input : Time series {s p (t j )}, p ∈ P , formatted as P = {x j } with x j = (s p (t j−d+1 ), · · · , s p (t j )) and y j = s p (t j+h ) d ∈ N (specifies sampling horizon), h ∈ N (specifies prediction horizon) c ∈ (0, 100) (percentage of data used for training) n o , n e , n r ∈ N, α o , α e , α r ∈ (0, 1] (parameters for function approximation) Let C contain c% of patients from P (drawn according to uniform prob. distr.) Set M = max{x j,k : x j ∈ C} and m = min{x j,k : x j ∈ C} for x j ∈ P do Rewrite x j = 2xj −(M +m) 2(M −m) end for ∈ {o, e, r} do for x j ∈ Q do Compute f (x j ) =F n ,α (C , x j ) end end this approach cannot be used for prediction based on data that was not included in the original data set we worked with. Indeed, this is one of the main motivations for our current method, which is not dependent on this requirement. Nevertheless, we include a comparison of the diffusion geometry approach with our current method for greater justification of the accuracy of our current method. Again, we used 50% training data and averaged the experiment over 100 trials, for prediction windows of 30 minutes, 60 minutes and 90 minutes and data sets D and J.
(iii) The Legendre polynomial prediction method followed in our 2013 paper [6]. In this context, the main mathematical problem can be summarized as that of estimating the derivative of a function at the end point of an interval, based on measurements of the function in the past. An important difference between our current method and the Legendre polynoimal approach is that with the latter, learning is based on the data for each patient separately (each patient forms a data set in itself, as explained in Section 1). Therefore, the method is philosophically different from ours. We nevertheless include a comparison for prediction windows of 30 minutes, 60 minutes and 90 minutes and data sets D and J.
(iv) The Fully Adaptive Regularized Learning (FARL) approach in the 2012 paper [10], which uses meta-learning to choose the kernel and regularization parameters adaptively (briefly described in Section 2). As explained in [10], a strong suit of FARL is that it is expected to be portable from patient to patient without any readjustment. Therefore, the results reported below were obtained by implementing the code and parameters used in [10] directly on data sets D and J for prediction windows of 30 minutes, 60 minutes and 90 minutes, with no adjustments for the kernel and regularization parameters. As explained earlier, this method is also not directly compared to our theoretically well founded method.
The percentage accurate predictions and predictions with benign consequences in all three BG ranges, as obtained using all five methods, for a 30 minute, 60 minute and 90 minute prediction window, are displayed visually in Figures 1-2. Tables 3-4 Figures 1-2, it is clear that our method far outperforms the competitors in the hypoglycemic range, except for the 30 minute prediction on data set D, where the diffusion geometry approach is slightly more accurate -but it is important to remember that the diffusion geometry approach has all the data points available from the outset. The starkest difference is between our method and Matlab's Deep Learning Toolbox -the latter achieves less than 1.5% accurate and benign consequence predictions regardless of the size of the prediction window or data set. There is also a marked difference between our method and the Legendre polynomial approach, especially for longer prediction horizons. Our method achieves comparable or higher accuracy in the hyperglycemic range, except for the 60 and 90 minute predictions on data set J, where the diffusion geometry and FARL approaches achieve more accurate predictions. The methods display comparable accuracy in the euglycemic range, except for the 60 and 90 minute predictions on dataset J, where Matlab's toolbox, the diffusion geometry approach and the FARL method achieve higher accuracy.

and the bar heights in
Tables 5-6 display the results when testing those methods that are dependent on training data selection (that is, our current method, the Deep Learning Toolbox and the diffusion geometry approach in our 2017 paper) on different training set sizes (namely, 30%, 50%, 70% and 90%). For these experiments, we used a fixed prediction window of 30 minutes, and training data was drawn randomly according to a uniform probability distribution in each case. While the Deep Learning Toolbox and diffusion geometry approach sometimes perform slightly better than our method in the euglycemic and hypoglycemic ranges, respectively, our method consistently outperforms both of these in the other two blood glucose ranges. The results are displayed visually as well in Figures 3-4.
As mentioned in Section 1, a feature of our method is that it does not require any information about the data manifold other than its dimension. So, in principle, one could "train" on one data set and apply the resulting model to another data set taken from a different manifold. We illustrate this by constructing our model based on one of the data sets D or J and testing it on the other data set. Building on this idea, we also demonstrate the accuracy of our method as compared to MATLAB's Deep Learning Toolbox when trained on data set D and applied to perform prediction on data set J, and vice versa, which is the only other method that is directly comparable with the philosophy behind our method. These results are reported in Tables 7 and 8. It is clear that we are successful in our goal of training on one data set and predicting on a different data set. The percentage accurate predictions and predictions with benign consequences when training and testing on different datasets are comparable to the cases when we are training and testing on the same dataset; in fact, we obtain even better results in the hypoglycemic  BG range when training and testing across different datasets.
Figures 5 and 6 display box plots for the 100 trials for the methods dependent on training data selection (that is, our current method, the Deep Learning Toolbox and the 2017 diffusion geometry approach) and prediction windows using data sets D and J, respectively. Our method performs fairly consistently over different trials.

Conclusions
In this paper, we demonstrate a direct method to approximate functions on unknown manifolds [7] in the context of BG prediction for diabetes management. The results are evaluated using the state-of-the-art PRED-EGA grid [14]. Unlike classical manifold learning, our method yields a model for the BG prediction on data not seen during training, even for test data which is drawn from a different distribution from the training data. Our method outperforms other methods in the literature in 30-minute, 60-minute and 90-minute predictions in the hypoglycemic and hyperglycemic BG ranges, and is especially useful for BG prediction for patients that are not included in the training set.     Table 8 Average PRED-EGA scores (in percent) for different prediction horizons on dataset J, with training data selected from dataset D.   We write ψ k (x) := h k (x) exp(−x 2 /2), x ∈ R, k ∈ Z + .