There is a growing trend of applying machine learning techniques in time series prediction tasks. In the meanwhile, the classic autoregression models has been widely used in time series prediction for decades. In this paper, experiments are conducted to compare the performances of multiple popular machine learning algorithms including two major types of deep learning approaches, with the classic autoregression with exogenous inputs (ARX) model on this particular Blood Glucose Level Prediction (BGLP) Challenge. We tried two types of methods to perform multi-step prediction: recursive method and direct method. The recursive method needs future input feature information. The results show there is no significant difference between the machine learning models and the classic ARX model. In fact, the ARX model achieved the lowest average Root Mean Square Error (RMSE) across subjects in the test data when recursive method was used for multi-step prediction.
The Blood Glucose Level Prediction Challenge with OhioT1DM Dataset [Marling and Bunescu, 2018] is aiming at developing powerful predictive models to predict the blood glucose level in 30 minutes. The dataset includes data collected from 6 anonymous patients with ID 559, 563, 570, 575, 588 and 591. For each patient, the following data were collected: a blood glucose level from continuous glucose monitor (CGM) every 5 minutes; periodic finger sticks blood glucose levels from; insulin doses, both bolus and basal; self-reported meal times with carbohydrate estimates; selfreported times of sleep, work, and exercise; and 5-minute aggregations of heart rate, galvanic skin response (GSR), skin temperature, air temperature, and step count.
This paper not only describes the detailed procedures of how the challenge results are obtained, but also conducts experiments to compare the performances of a wide spectrum of machine learning algorithms. In the light of No Free Lunch Theorem [Wolpert, 1996] 1, it is worth trying out various learning algorithms on a given problem. In the meanwhile, plenty of literature have reported the blood glucose level prediction results separately with different datasets using either classical time series approaches [Gani et al., 2009; Eren-Oruklu et al., 2009; Sparacino et al., 2007; Turksoy et al., 2013; Wang et al., 2014; Xie and Wang, 2017] and machine learning approaches [Zecchin et al., 2012; Plis et al., 2014; Mirshekarian et al., 2017; Mhaskar et al., 2017; Fox et al., 2018]. However, it is hard to compare the performance of different approaches since different datasets were used by different papers. Even with the same dataset, different features or regressors might be used to train the model, which profounds the comparison of different algorithms. We hope the comparison results listed in this paper is able to provide some empirical insights of algorithm selection on this particular blood glucose level prediction problem. The algorithms that are implemented for comparison include the classic autoregression with exogenous inputs (ARX) model, Huber Regression, Ridge Regression, Elastic Net Regression, Lasso Regression, Support Vector Regression with Linear Kernel, Support Vector with Radial Basis Kernel, Random Forest, Gradient Boosting Trees. In addition, the results of two deep learning models, a vanilla Long-Short-TermMemory (LSTM) Network and a Temporal Convolution Network (TCN) [Bai et al., 2018], are also reported. It is worth pointing out the LSTM Network model has its special mechanism to handle the memory of the sequence, which diverges from the typical autoregression manner used by the other approaches. Keep in mind that the results of LSTM Network might not be a fair comparison with the other models. 2
Here we introduce the general autoregression model structure which is suitable for most of the machine learning algorithms (except Recurrent Neural Network which maintains hidden states to keep memories) on time series prediction problems. In the end of the section, we introduce two methods to make multi-step predictions. To make the results of this paper comparable with other results, we will report the results using both multi-step prediction methods.
1There is no such a learning algorithm that outperforms any other learning algorithms in every problem set. (1) The general autoregression model relies on the histories of both the target values and input feature values to make kstep-ahead predictions. The model structure is represented as follows.
y(t + k) =f (y(t); :::; y(t
p + 1); X(t); :::; X(t q + 1)) + e(t) where y(t) is the target time series to predict, X(t) is a multivariate time series containing the input features used for prediction. p and q are the regression orders for the target sequence and the input sequence. The noise term e(t) is assumed to be identically independently distributed. Now the time series prediction problem is formulated as a general regression problem. Function f ( ) is the actual model that varies among different machine learning algorithms. All the non-deep-learning models implemented in this paper fall into this category. By fixing the regression orders p; q and varying the model type of f ( ), fair comparison of different machine learning algorithms can be made. Note that for one-stepahead prediction problem, when the ordinary least squares (OLS) algorithm is used to fit f ( ), the predictive model becomes an ARX model with input order q and output order p. 2.2
Recurrent Neural Network (RNN) and Convolution Neural Network (CNN) are the two major types of deep neural networks that achieved notable improvements recently in multiple machine learning tasks. We selected a vanilla LongShort-Term-Memory (LSTM) RNN structure [Hochreiter and Schmidhuber, 1997] and a Temporal Convolution Network (TCN) [Bai et al., 2018] as representatives to benchmark their performances in this Blood Glucose Level Prediction problem.
The vanilla LSTM network implemented in the paper is simply cascading multiple layers of LSTM cells with a fullyconnected linear layer in the end to produce the predicted blood glucose value k step ahead. Unlike the general autoregression model, the prediction of the general RNN depends on both the current measurements and the recurrent hidden state: y(t + k) = f (h(t); y(t); X(t)) h(t + 1) = g(h(t); y(t); X(t)) (2) RNN relies on the hidden state to keep its memory of the sequence. This fundamental structural difference between the general autoregression model and RNN model makes it difficult to make fair comparisons of their performances.
The Temporal Convolution Network (TCN) we implemented is identical to the network introduced in [Bai et al., 2018]. TCN adopted the key structures from state of the art networks WaveNet [Van Den Oord et al., 2016] and ResNet [He et al., 2016]. It utilizes the dilated causal convolution of WaveNet to obtain long term memories of the input sequence, and uses the residual connection inspired by ResNet to make training deep network easier. Similar to the vanilla LSTM Network, a fully-connected linear layer is appended in the end to make the final prediction.
Note that TCN is stateless when making predictions. The target prediction at time t + k is still a function of its receptive field, i.e.
y(t + k) =fT CN (y(t); :::; y(t r + 1); X(t); :::; X(t r + 1)) + e(t) (3) where r is the receptive field length of the TCN model. For a n layer TCN with convolution kernel size k and dilation factor d, the receptive field r is equal to (k 1)dn [Bai et al., 2018]. Hence, TCN is still under the framework of the general autoregression model, which makes fair comparison possible with other machine learning algorithms. 2.3
There are two major strategies to obtain multi-step prediction: recursive method and direct method. To apply the recursive method, first a model is trained to predict the target value one step ahead. Then the predicted one-step ahead value is reused to produce the make the prediction of the next step. By iterating this procedure, one is able to obtain multi-step prediction. It should be noted that predicted input feature values X^ (t+k) is also needed to make recursive predictions. To simplify the prediction procedure, the ground truth values of input features values X(t + k) was used instead of the predicted values in the paper, so the prediction error only reflects the uncertainties that the model fails to capture (if predicted values X^ are used, the final prediction error is the uncertainty of X plus the model uncertainty).
Direct method treats the k-step-ahead value y(t + k) as target values directly and only uses the input feature values till current time t. There is no future information leakage using this method. The results using direct method is not comparable with the results using recursive method in this paper. 3
The data in the OhioT1DM Dataset is highly unsynchronized, since the data are collected from multiple devices, and some of the data is logged by patient manually. In addition, missing data exists in the data collected from the CGM sensors and the Basis Peak fitness band. Therefore, the data cleaning procedures can be split into two stages: resample stage to align all the data onto the same time grid, and imputation stage to fill the missing data.
At the resample stage, a time grid with 5 minute sample period was derived based on the start and end timestamps of the CGM blood glucose signal. Then all the signal were resampled to the time grid in a “forward filling” manner, which means if the value is missing at some time point on the designated time grid, then the closest value before the time point will be used. In addition to the “foward filling” rule, the same value can only be used once at most in the resampling process. As long as the closest previous value has been used once for resampling, the following missing values will just be tagged as missing. By resampling the data in this way, we are able to reindex the signals onto the same timestamps and preserve the locations of the missing data at the same time. To perform this cleaning process, reindex function in the pandas package of python was used with method set to ffill, and limit set to 1. To be noted that, this resampling strategy was only applied to CGM sensor data and data collected from the Basis Peak fitness band. For the other event based data like basal bolus insulin, meal and exericse, the resample strategy is based on the context provided by the data description [Marling and Bunescu, 2018]. For example, the <basal> and <temp basal> data only contains basal rate change time and the changed value. We combined the information provided by both to obtain the continously sampled insulin basal rate. Also notice that insulin bolus is in unit of dose in the original data. It was converted to rate (dose/min) by dividing “normal” or “square normal” type bolus by 5 minutes (sample period) and dividing “square dual” type bolus by its duration during the resampling process. Figure 1 shows the resampled data of subject 559. At the imputation stage, each signal with missing data is imputed by a Kalman Smoothing technique described in [Moritz and Bartz-Beielstein, 2017]. The kernel model used by the Kalman Smoothing technique is a structural time series model introduced in [Commandeur et al., 2011]. It is important to keep in mind that this Kalman Smoothing imputation technique requires the future information to perform the imputation. It cannot be used to impute data in real time. In practice, na.kalman function in the ImputeTS package of R was used with method set to “StructTS”. Figure 2 shows the CGM blood glucose level before and after imputation. Same cleaning procedures were performed in the final testing dataset as well. 4
engineering we made was adding the insulin basal rate and bolus rate (the raw bolus dosage data was converted to dosage rate during the data cleaning stage) to get the total insulin delivery rate2. To benchmark different learning algorithms, same input features X(t) were used for training, which are the derived total insulin delivery rate, meal sizes, and the heart rate measured by the Basis Peak fitness band. 5
While taking the data challenge, only the training data is available for model development. To evaluate the performances of each learning algorithm, the cleaned training data was split into sub-training data and sub-testing data. We only used the sub-training data to train the models. The sub-testing data was leveraged to pick the model for final submission of the challenge. The model that achieved the best performance on the sub-testing dataset were used to make the final prediction of the testing dataset of the challenge. In our case, the Linear Regression (ARX) model was selected for final challenge result submission.
For subject 559, 563, 570 and 588, the sub-training data takes the first 80% of the original data, while the rest is treated as the sub-testing data. For subject 575 and 591, noticeable amount of continuously missing data is observed between date Dec 26 to Dec 28 and Dec 26 to Jan 5 respectively. Therefore, the dates of missing data are used as split point for subject 575 and 591. 6
For each learning algorithm, we trained two models for each subject, one model for 1-step-ahead (5-minute-ahead) prediction and the other for 6-step-ahead (30-minute-ahead) prediction. The 1-step-ahead predictive model was used in the recursive method to produce 6-step-ahead predictions.
For training non-deep-learning models, the sub-training data is further split into 4 folds. Then three cross-validation datasets were generated: 1st fold for training and rest for Feature engineering is not a focus of this paper. Hence little effort was spent to search for useful features. The only feature 2The units of basal and bolus insulin are assumed to be the same, since they are not provided by the dataset. Huber
Ridge ElasticNet
Lasso SVR-LinearKernel SVR-RadialBasis
RandomForest GradientBoostingTrees
Hyperparameters validation; 1st and 2nd folds for training and rest for validation; first 3 folds for training and the last 1 fold for validation. Grid search was performed for each learning algorithm to fine-tune the hyperparameters. The model with the lowest average mean square error across all three cross-validation datasets was saved, and was used later on to generate multistep predictions (6 steps to get 30-minute-ahead predictions in this case) on the sub-testing dataset and the final challenge testing dataset. The hyperparameters used in each model is listed in Table 1 and Table 2. Before fitting the model, the blood glucose, insulin, meal and heart rate data was normalized to the scale of 0 to 1. Note that the normalization only utilized the information of data to fit, which prevented information leakage when producing the validation scores. Exactly the same normalization operation was performed on the sub-testing data and the final challenge testing data. The regression orders p and q were set to 12 (1 hour) for all the learning algorithms.
While training the two deep learning models, the sub-training data is further split into 5 folds with the first 4 folds for training and the last fold for validation. The validation data was used for early stopping. In practice, the training algorithm checks the validation score at the end of each epoch, and only saves the model when it achieves the best validation score so far. The training and validation data was also normalized to the scale of 0 to 1 in the same way of training non-deeplearning models.
At the training stage, the original long sequence data of each subject was truncated in a sliding window fashion so that each subject will have multiple training examples. Here we set the sliding window size to 576 samples (48 hours), and step size was 12 samples (1 hour). Adam optimization with learning rate 0.001 was used in training both models. The coefficients used for computing running averages of gradient and its square were set to 0.9 and 0.999 respectively. Learning rate is shrunk to its 1=10 if the model has been trained for more than 10 epochs and the validation loss is larger than the loss in the last 3 epochs. The gradient norm was clipped at 1. Batch size was set to 1. All the deep learning networks were trained for 20 epochs.
The vanilla LSTM Network has 3 hidden LSTM layers with hidden size of 50. The TCN model has 2 layers of TCN blocks (referred to TCN2 in Table 4 and Table 5) with kernel size of 4 and dilation factor of 2, so that it has the same regression order of 12. To investigate the performance of the TCN model with a deeper structure, another TCN model with 10 layers of TCN blocks (other parameters were the same as the shallower one) was implemented (referred to TCN10 in Table 4 and Table 5). Dropout of 0.5 was applied after each TCN layer. It should be noted that both the vanilla LSTM Network and the 10-layer TCN model might not be a fair comparison with other models. 7
The final challenge testing dataset was first cleaned under the exactly the same procedures as the training data. With the cleaned challenge data, two sets of 30-minute-ahead prediction were produced by the 1-step-ahead models and the 6-step models using recursive method and direct method respectively. Table 3 shows the test points change after each data cleaning step, and the final test points used to produce the challenge results. Though the target blood glucose signal was imputed during the cleaning stage, the imputed values were ignored when computing the final RMSE scores. Since the regression order p = q = 12, and the model is making 6-step-ahead prediction, the first 12 + 6 1 = 17 test points were ignored as well when generating the final results. That is why the final test points used are 17 points less than the test points after resampling. It should be noted that 1 test point was removed from the data of subject 570, and 8 new test points were introduced to the data of subject 575 during the resampling stage. The results of 575 might not be comparable in the challenge.
Subject Original Resampled Imputed Test point used 570 The final test scores for each model using recursive method and direct method are reported in Table 4 and Table 5. The
Linear Regression (ARX) ElasticNet GradientBoostingTrees Huber Lasso RandomForest Ridge SVR-LinearKernel SVR-RadalBasisKernel VanillaLSTM* TCN2 TCN10*
Linear Regression ElasticNet GradientBoostingTrees Huber Lasso RandomForest Ridge SVR-Linear Kernel SVR-RadalBasisKernel VanillaLSTM* TCN2 TCN10* 23.20 18.36 19.28 22.60 19.09 19.89 21.35 18.36 18.69 19.13 69.57 19.86 20.74
559 23.20 18.30 18.71 20.53 18.47 18.86 21.07 18.30 18.25 18.19 21.23 19.86 19.15 20.71 19.02 18.51 19.74 20.39 18.85 20.04 19.02 19.71 20.08 56.30 18.91 19.50
563 20.71 19.11 18.51 19.79 19.08 18.51 20.22 19.11 19.29 19.12 18.62 18.91 21.24 19.05 16.03 16.83 17.61 16.16 17.59 16.90 16.03 15.73 15.84 46.28 19.60 17.67
570 19.05 16.02 18.95 18.38 16.10 18.95 18.05 16.02 15.69 15.67 17.12 19.60 17.60 25.63 23.90 24.08 23.83 25.95 24.38 24.60 23.90 25.00 25.59 90.47 23.46 27.55
575 25.63 23.85 23.92 23.80 25.24 23.97 24.47 23.84 24.61 24.61 25.09 23.46 24.39 21.99 18.25 19.13 18.82 18.65 20.17 19.72 18.25 17.90 17.66 71.72 20.72 20.74
588 21.99 18.20 18.66 19.32 17.98 18.91 19.71 18.20 17.82 17.49 19.81 20.72 20.21 24.64 21.99 22.04 22.30 24.48 22.25 22.31 21.99 22.80 22.23 73.09 25.35 23.33
591 24.64 22.25 22.12 22.53 22.69 22.12 22.76 22.25 22.65 22.12 22.47 25.35 24.65
Mean 2.25 2.85 2.62 2.45 3.72 2.44 2.61 2.85 3.38 3.45 15.20 2.53 3.46 SD 2.25 2.89 2.30 2.06 3.38 2.27 2.29 2.89 3.31 3.27 2.85 2.53 2.84 zero order hold approach (use current measurement as the prediction in 30 minutes) is used as a baseline performance.
The results using recursive method indicate that the classical ARX model performs the best in average. Ridge Regression achieved almost the same results as ARX because the coefficient of regularization penalty was set to 0.0001 as a result of grid search, which makes the Ridge Regression nearly identical to the Linear Regression used by ARX. Nonlinear models like SVR with Radial Basis Kernel, Random Forest, Gradient Boosting Trees and the two types of Deep Neural Network did not outperform other linear models. Performance degraded for some subjects when using a deeper 10-layer-TCN structure compared to the 2-layer-TCN.
It is noticed that the RMSE of the vanilla LSTM Network using recursive method is significantly higher than other models. Further investigation found that the one-step-ahead prediction of the vanilla LSTM is still accurate, while the recursive multi-step prediction accumulated undesirable errors resulting in large oscillation of the prediction. The vanilla LSTM Network potentially overfits the training data. The one-step-ahead performance of the vanilla LSTM Network is listed in Table 6.
When using direct method to make multi-step predictions, SVR with Radial Basis Kernel had the best average performance, while the Linear Regression is still very competitive with the second best average performance. Again Ridge Regression selected 0.0001 penalty on the regularization terms during the grid search, so its performance was nearly identical to Linear Regression. Unlike the recursive method, the direct method does not prone to error accumulating issues, and the vanilla LSTM Network is doing a decent job in the multi-step prediction compared to its poor performance in recursive method. It is worth pointing out that TCN models failed to outperform the baseline model for some subjects. 9
This paper benchmarked multiple popular machine learning algorithms including the classical ARX method and 2 modern deep neural networks on the OhioT1DM Dataset. The results indicate that the classical ARX model achieved the lowest RMSE in the final testing dataset when making multistep prediction in a recursive manner. The Linear Regression is also a good choice when making predictions with the direct method, since there is no need of hyperparameter tuning in Linear Regression and short training time can be achieved. Though the predictive models achieved decent performance in glucose level prediction, whether the predictive model captured the correct relationships between blood glucose and the input features like insulin, meal and heart rate is the still an important aspect to investigate in the future.