A machine learning-based method for blood glucose level prediction thirty and sixty minutes in advance based on highly multiclass classification (as opposed to the more traditional regression approach) is proposed. An advantage of this approach is the possibility of modelling and visualising the uncertainty of a prediction across the entire range of blood glucose levels without parametric assumptions such as normality. To demonstrate the approach, a long-short term memory-based neural network classifier is used in conjunction with a blood glucose-specific data preprocessing technique (risk domain transform) to train a set of models and generate predictions for the 2018 and 2020 Blood Glucose Level Prediction Competition datasets. Numeric accuracy results are reported along with examples of the uncertainty visualisation possible using this technique.
Maintaining blood glucose level (BGL) in the normoglycemic range
is a significant challenge for patients with type 1 diabetes (T1D).
Traditionally, patient BGL self-management is achieved using finger
stick blood samples, testing strips and glucose meters (see [
In this paper, the problem of forecasting BGL thirty and sixty
minutes in advance is considered using the 2020 BGL Prediction
Challenge [
Here, an alternative approach is taken: instead of treating BGL
forecasting as a regression problem, it is instead viewed as a
classification problem. This is achieved by dividing the range of
possible BGL values into 100 bins equally spaced in the risk domain
[
In order to transform the problem of BGL forecasting into a classification problem, a method of breaking the BGL range into sensible classes is required. This is not trivial because the range of BGL values is continuous, and the sizes of clinically-relevant subranges varies non-linearly. For example, a small change in the hypoglycemic part of the BGL range may be highly significant clinically but an equivalent change in the hyperglycemic part of the range may be considered insignificant.
One option is to use the five ranges proposed by Danne et al. [
Because of the inherent non-linearity of the BGL range, an
approach called the risk domain transformation, first proposed by
Kovatchev et al. [
As can be observed in the figure, risk domain values typically spans a range from approximately -2 to just over 2, and most normoglycemic readings lie more or less in the range range [ 0:9; 0:9]. A risk value of 0.0 corresponds to the BGL of 112.5 mg/dl, which is considered the point of least risk. An advantage of the risk domain is that the hypo- and hyperglycemic ranges now have equal size and significance, which reduces the chance of bias in statistical analysis (e.g. due to larger absolute error sizes in the hyperglycemic range). r(xt) = yt = 1:509 log(xt)1:084 5:381 r 1(yt) = xt = exp
yt 1:509 + 5:381 For dividing the BGL range into classes, the following exact procedure is used. The risk domain range [ 2; 2] is considered and 100 bin midpoints are placed on it. The bin midpoints are denoted y1 ; y3 ; y3 : : : y100 with y1 = 2 and y100 = 2. The remaining bin midpoints are equally spaced along the risk range between y1 and y100. This ensures that for the smaller hypoglycemic range, the bin sizes will be scaled properly in size and that there will be a proportionate number of bins (and therefore classes) in each subrange of the BGL scale.
Next, in order to a assign a new CGM reading x (in mg/dl) to a bin, its corresponding risk value y = r(x) is computed using Equation 1. The reading is then assigned the bin with the closest midpoint. This means that the split points between bins do not need to be calculated explicitly, and if a reading is outside the risk range (either x < 2 or x > 2) then it will be assigned to one of the bins at the ends of range. However, this tends not happen often since the significant majority of readings in the competition datasets lie on the [ 2; 2] range.
Finally, once its class is determined, the reading x is transformed into a one-hot encoded vector (0; 0 : : : 0; 1; 0 : : : 0; 0) of length 100 where the single 1 in the vector corresponds to the bin that x is assigned to.
To summarise the process, the predictors X for the model comprise a time series of risk-transformed CGM readings, and the target Y is a one-hot encoded class vector of dimension 100. Predictions are therefore numeric vectors, e.g. (0; 0 : : : 0:25; 0:6; 0:3 : : : 0; 0). Note that whichever type of model is used, the values should be positive and sum to unity so that they can be interpreted as normal probabilities.
To evaluate this idea, 100-class classification experiments using a neural network as a predictive model were performed. Figure 2 depicts the particular neural network architecture used here.
LSTM(units=12)
Dense(units=50, act=RELU)
Dense(units=100, act=Softmax)
The inputs to the neural network are twelve risk value readings, representing the past sixty minutes of BGL variation as sensed by the CGM (assuming that readings occur every five minutes). Since this is a time series, a long-short term memory (LTSM) layer with twelve units is used for initial input processing. Next, the LSTM output is flattened and passed through two dense fully connected layers for further processing. Both of these layers apply batch normalisation first, which ensures faster convergence times and stability during training. Finally, the last dense layer predicts the class and applies a softmax activation function to ensure that the output is a probability vector.
To train each instance (one per patient) of this neural network, the
ADAM optimiser [
To generate data for training the neural network, a strict approach was taken towards missing data and all examples with gaps or time discrepancies (e.g. readings not exactly five minutes apart) exist were excluded. Therefore, to generate one example for the thirty minute (t+30) forecasting problem, it is required that nine consecutive readings (from the start of the example up to and including the prediction target) exist, and for the sixty minute problem (t + 60), twelve consecutive readings were required. No missing value imputation was performed. As a result, the number of test examples varies slightly depending on the forecasting horizon: i.e. the number of 2020 dataset test examples is 2,743, 2,579, 2,177, 2,185, 2,393 and 2,624 for the 30-minute horizon and 2,689, 2,531, 2,111, 2,113, 2,297 and 2,582 for the 60 minute horizon respectively.
With the neural network architecture and training data construction approach described, the final aspect of methodology to be described is the way that numeric point predictions were generated for the competition purposes (which require point predictions). A simplistic approach is allow only bin midpoints (i.e. y1 : : : y100) to be predictions, and select the bin/class with the highest probability. Initial tests showed that this technique had low accuracy. Instead, a more sophisticated approach is to calculate the expected BGL value, as described by the following equation: x^t+n = f (mn; yt 55; yt 50 : : : yt 5; yt)
100
= X p(yi jmn; yt 55 : : : yt)r 1(yi )
i=1
(3)
where n 2 f30; 60g, mn is the neural network for the current patient
with forecasting horizon n, yt is the risk-transformed BGL level at
time t (yt = r(xt)), y1 : : : y100 are risk bin midpoints, f () represents
the application of mn to the observed CGM values, and x^t+n is the
expected value or prediction at time t + n. Since the probabilities
across the 100 bins sum to 1, the resulting point prediction will be
scaled correctly. Initial experiments showed that this expected value
approach produced accurate estimates. Source code is available [
Two rounds of experiments were performed. In the first round, neural networks models were independently trained and tested for each of the twelve patients from both the 2018 and 2020 competition datasets. There was no sharing of information between models. In the second round, models for the 2020 patients only were trained, and the training data for each patient included all of the training data from the 2018 competition in addition to the specific 2020 patient’s training data. An individual patient’s own training data was therefore a small subset of his/her full training dataset. To account for this imbalance in the second round, samples from the target patient were re-weighted by a factor of six compared to the sample weights from the other patients.
Results for the first round of experiments are given in Tables 1 and 2. The first table is mean absolute error (MAE) results, and the second table gives root mean squared error (RMSE) results. Units are mg/dl. Using both metrics, predictive performance is comparable between the two datasets, with patient 540 from the 2020 dataset being the most “difficult” patient to predict. Conversely, the patient with lowest forecast error is patient 570 from the 2018 dataset.
Results for the second round of experiments are given in Tables 3 and 4. It can be observed that the additional training data leads to a very slight improvement in accuracy. Performing a paired t-test across the six 2020 patients reveals that the improvement in MAE is significant even albeit small (average improvement for thirty minute forecasting is 0.45 with significance p = 0:000056, and for sixty minute forecasting it is 0.82 with significance p = 0:008096).
More interesting are the prediction plots that can be generated when making a forecast using the classification approach. Figures 3-6 depict some probability densities produced by the model when making four different predictions. For each figure, the point prediction generated using the expected value computation (Equation 3) is shown as red line.
The tidiest example is Figure 3, which depicts a single-peaked distribution with the forecast coinciding with the peak of the distribution. While this class of forecast is common, it is not the only type of distribution that is output from the model.
Figure 4 shows a two-peaked distribution with an expected value between the peaks. While the expected value is a useful point prediction, the dual peaks are also useful information since it can be clearly observed that the two most likely outcomes are normoglycemia vs. a state most likely in level I hyperglycemica, although the model is not certain.
Figure 5 shows a skewed distribution a significant mass of the probability distribution is at the upper end of the range, but the expected value is closer to the middle of the range (albeit still in the hyperglycemic range). In this case, the expected value underestimates the true risk to the patient at the current time. Again, this is clearly noticeable.
Finally, Figure 6 depicts a prediction with considerable noise in the distribution. While the forecast is normoglycemic, there is significant mass in the hypoglycemic range. Therefore it could be concluded that although the point prediction is reasonably, there is still significant risk of hypoglycemia.
Further analysis was performed with respect to the variance of the prediction distributions. It was found that for most patients, the distribution of variances is skewed to the left indicating that on average most predictions are more certain (more like Figures 3 and 4) than uncertain. However, more analysis needs to be done on this point. 4
This paper describes a system for forecasting BGL at thirty and sixty minutes in advance. This main distinctiveness of this approach is the adoption of a highly multi-class classification-based technique and use of a domain-specific transform for normalising BGL values (opposed to more traditional min/max scaling or standardisation). The ability to visualise non-parametric probability distributions accompanying predictions as a meaningful context is a clear advantage.
To test the proposed method with real patients, we will use a
system known as SmartCGMS [
Our specific approach will be to transform the best trained keras/tensorflow-based neural network into a hard-coded and constant feed-forward neural network in C++. This will enable efficient deployment and computation on low-power devices such as insulinpump controllers, while we can still train the original neural network using high-performance computers. As a result, a flow in which a neural network is continuously learned from patient BGL measurements, providing personalised BGL predictions, can be established.
This publication was partially supported by the project LO1506 of the Czech Ministry of Education, Youth and Sports.