Recent developments in diabetes care technologies, like wearable devices for continuous glucose monitoring, have made it possible for patients to access relevant data and properly treat diabetes mellitus diseases. Thanks to machine learning, moreover, it has been possible to predict the trend of the level of glucose in the blood in real-time, thus preventing diseases. In this research, we propose to compare the two main state-of-the-art algorithms for time series forecasting, namely Autoregressive Integrated Moving Average (ARIMA) models and Recurrent Neural Network (RNN) models. In particular, we propose an auto adaptive algorithm for ARIMA models' parameters selection based on the Augmented Dickey-Fuller test for stationarity and autocorrelation and partial autocorrelation functions for autoregressive and moving average orders and we show that somehow this would make the ARIMA model preferable to RNNbased models in an online learning scenario for time series forecasting, with a root mean square error of 1.11 mmol/l. We also propose some considerations related to a Google Cloud based infrastructure to host an online learning application, comparing the performance of the tested models also with respect to their training time and their scalability and maintenance, showing that ARIMA models clearly require a lighter infrastructure and a less complex pipeline for managing model life cycle with respect to RNN-based models.
Diabetes mellitus is a metabolic disorder that causes blood glucose level to deviate from
normal values and can lead serious health complications and in some cases, if not
properly treated, even death [
In descriptive statistics, a time series is defined as a collection of constant points in time
and expresses the dynamics of a certain phenomenon over time. The analysis of time
series has two main practical applications: on the one hand, to provide an
interpretation of a phenomenon, identifying components of trends, cyclicality and/or seasonality;
on the other, to predict the future trend of the phenomenon [
The goal of this research is that of defining an online learning algorithm for glycemic index forecasting. The research is, indeed, driven by business needs since the developed algorithm is intended to be used within a system that aims at monitoring some vital parameters of bus drivers during their working hours in order to prevent diseases. In particular, data acquisition takes place thanks to different devices that monitor (with different frequencies) the blood glucose level in terms of glycemic index, the lipid profile and HbA1c, the blood pressure and heart rate. Some of these variables are acquired manually by a specific device once or twice a day. For these variables a prediction is not needed, but the collected values are visualized in order to monitor their trend. As far as the glycemic index is concerned, indeed, it is collected by a wearable device called Glunovo CGM every 3 minutes, and for this variable the goal is to predict the trend of the series in the near future, in order to alert the driver in case the trend is going out of the predefined safe boundaries. Thus, the final system must include: a component for the prediction of the glycemic index signal in the near future; an alerting system capable of reporting dangerous situations; a dashboard that can be consulted by the doctor who, based on the predictions made and the other values collected, decides whether to alert the driver or not. In this research we implement an algorithm that, at each new value registered by the Glunovo CGM device, analyzes the historical series up to that instant and provides a real-time forecast of the trend of the glycemic index for the following period with an horizon of 30 to 60 minutes. Moreover, the algorithm must be able to adapt to new registered values in an online learning fashion and it must be suitable for a huge number of drivers. 3.1
In order to design the online learning framework, we use an open source dataset called
D1NAMO suitable for the analyzed case. The dataset is described by F. Dubosson et
al. in [
yt = t + t (1)
A series is called autoregressive of order p, AR(p), if the value of the series at instant t is a linear combination of the values of the signal from instants t 1 to tp.
yt = 1yt 1 + 2yt 2 + ::: + pyt p + t (2)
A series is said to have a moving average of order q, M A(q), if the value of the series at instant t is a linear combination of the values that describe the trend of the series from instants t 1 to tq .
A series is called autoregressive moving average of orders p, q, ARM A(p; q) if it is both autoregressive of order p and moving average of order q.
yt = 1yt 1 + 2yt 2 + ::: + pyt p + t + 1 t 1 + 2 t 2 + ::: + q t q (4)
AR(p); M A(q) or ARM A(p; q) models can only model stationary time series. If the joint probability distribution function of fyt s; :::; yt; :::; yt+sg is independent of t for all s, a series yt is strongly stationary. Thus, for a strong stationary process the joint distribution of any possible set of random variables from the process is independent of time. However for practical applications, the assumption of strong stationarity is not always needed and so a somewhat weaker form is considered. A stochastic process is said to be weakly stationary of order k if the statistical moments of the process up to that order depend only on time differences and not upon the time of occurrences of the data being used to estimate the moments. A homogeneous non-stationary series can, however, be reduced to a stationary series by introducing an appropriate degree of differentiation d. Using this degree of differentiation d, if the series is also autoregressive with moving average, the series is called autoregressive integrated with moving average of order p; d; q ARI M A(p; d; q). It is important to note that neither strong nor weak stationarity implies the other. However, a weakly stationary process following normal distribution is also strongly stationary. To verify if a series is stationary it is possible to use some common statistical test like the one given by Dickey and Fuller. To determine the optimal parameters for an ARMA model, it is necessary to carry out the ACF and PACF analysis. These statistical measures reflect how the observations in a time series are related to each other. The ACF is used to measure the correlation between the current observation and an observation at lag k (using the Pearson coefficient of correlation by assuming that each variable is distributed as a Gaussian), while the PACF is used to measure the correlation between the current observation and an observation at lag k, after removing the effect of any correlation due to observations at intermediate lags (i.e. at lags < k ). For each value of the ACF and PACF also a confidence interval is computed to measure the significance of that computed correlation. For modeling and forecasting purposes it is often useful to plot the ACF and PACF against consecutive time lags. These plots help in determining the order of AR and MA terms. Some empirical considerations, indeed, relates to the following facts: if the ACF decreases slowly and the PACF decreases very quickly after p steps, the series is AR(p); if the ACF decreases very fast and the PACF decreases slowly after q step, the series is MA(q); if both the ACF and the PACF decrease slowly, the series is ARMA. By analyzing the graph shown in Figure 2, it is possible to see the ACF and PACF of a time series for 2 degrees of differentiation. Auto adaptive ARIMA parameters optimization algorithm Given the previous considerations, in this research we define an adaptive logic for optimizing the parameters of the model p, q and d so that, for each new data point, the best model is recalculated in real-time on the basis of the historical data in an online learning fashion. In this way, the trend of the series is predicted for each new data point. Since prediction must be taken in real-time, it is important to quickly compute the selection of the parameters. The algorithm we designed works iteratively and is reported in pseudo-code in the following. At each iteration it takes the historical series (values) and calculates the differential of the series starting from 0 by incrementally increasing the value of the differential to be calculated. Once the differential has been calculated, the ADF statistical test is performed. This test verifies the null hypothesis that there is a root in the time series. The alternative hypothesis is the stationarity hypothesis. Through this test, it is therefore possible to determine with statistical significance whether the series is stationary or not. The test is verified by considering the value of the p-value and the value of the calculated ADF statistic. If the p-value < 0.05 and the calculated statistic is less than the first critical value returned by the test, then the series is considered stationary, otherwise the series is considered non-stationary. Once the test has been performed, if the null hypothesis holds, the degree of differentiation of the series is increased and the test is repeated, if the alternative hypothesis holds, instead the series is considered stationary and the algorithm proceeds in the determination of the parameters p and q. To determine such parameters, we use the calculation of the ACF and PACF functions. By following the empirical considerations previously described, to detect p and q we compute the cut-off of these series and take the first value of the function which is outside the confidence interval (considering an alpha of 0.05). We therefore consider the maximum value that p can take, p max, as the cut-off value of the ACF function and the same for q, q max, as the cut-off value of the PACF function, i.e. as the number of detections of the ACF and PACF functions that are outside the 95% confidence interval. The choice to exclude the values that are outside the confidence interval is determined by the fact that those values do not have an effective statistical significance in determining the trend of the series. For this reason, the hypothesis is that the p and q parameters should not be higher than these cut-off values, however it is not guaranteed that the cut-off values are the optimal model values. Indeed, once the p max and q max parameters have been determined, to guarantee the optimal choice of the model, a greedy search of the parameters is carried out, keeping d fixed and making p vary from 1 to the ACF cut-off value and q from 0 to the cut-off value of the PACF. For each model trained, the AIC is calculated. This parameter is an estimator of the prediction error and, consequently, of the relative quality of the statistical model calculated on the input data. Indeed, given a collection of models for a dataset, the AIC estimates the quality of each model relative to the quality of the others. In this way, the AIC provides an estimate to make the selection of the models. In this way it is possible to select the best possible model and with each new data the model can be retrained and used for prediction.
Algorithm 1: ARIMA Optimal Parameter Selection
Input: The time series values ts val to train the ARIMA model Output: The optimal parameter p, d, q of the ARIMA model for the given time series d 0; check False; while check=False do diff val differentiate ts val with order d; p value, ADF stat, critical val compute ADF test over diff val; if p value<0.05 and ADF stat<critical val then acf val compute ACF over diff val; pacf val compute PACF over diff val; p max compute cut off of acf val; q max compute cut off of pacf val; check True; else end d
increment d by 1; end best p 0; best q 0; aic val 1; for all p~ 2 [1; p max] do for all q~ 2 [1; q max] do arima model compute ARIMA(p~, d, q~); aic compute AIC for the arima model; if aic < aic val then best p p~; best q q~; end end 4.2
An RNN is a class of artificial neural network in which the output values of a layer of
a higher level are used as input to a layer of a lower level [
st = f (U xt + W st 1 + bs) (5)
In the above equation, f is the activation function of the hidden layer and bs is the bias vector of the hidden unit. The output layer has P units ht = (h1; h2; :::; hP ) calculated as below, where V is the matrix of weights that connect the hidden layer with the output layer, g is the activation function of the output layer and bh is the bias vector of the output unit.
h = g(V st + bh) (6)
Thus, RNN represents a non-linear equation of state iterable over time. At any time,
the hidden state provides a prediction of the output based on the input vectors [
In order to define the best solution for the given problem, we decided to compare
such LSTM-based models trained on Google Cloud Vertex AI leveraging a
hyperparameter tuning tool based on Google Vizier known as Google HyperTune [
As mentioned in Chapter 3, a public dataset (D1NAMO) is used to train and test the models. For all tested models, two prediction time horizons are compared, one at 30 minutes and one at 60 minutes. For both cases the training is carried out considering the measurements in the previous 8 hours. Thus, the time series is divided into sub-serieses of 8 hours and divided into train and test according to the 80-20 split. Since in the dataset the detections have a sampling period of 5 minutes, the prediction is made for 6 and 12 subsequent points and the training is done using 96 samples for each series. To measure the results we use RMSE defined as:
T X(y~t t=1 yt)2jj (7)
As far as the ARIMA model is concerned, it is trained with the last 96 samples of training data and over these the prediction is carried out (to be compared over the first 6 or 12 values of the test set). Then for each subsequent point with respect to the last 96 the ARIMA model is re-trained and the prediction is compared with the test data. This is done for 7 patient-related time series out of the 9 available, given the relative scarcity for 2 of these 9, and till all the test data are covered. Results are then averaged. As far as LSTM is concerned, we train a model for each patient over the training data, and each model is used to compute the prediction for subsequent values of the time series in the test set. Also in this case the results are averaged. Glucose values recorded in the reference dataset are measured in mmol/l, which has a conversion rate in mg/dl s of 18.02. Furthermore, given the need to make predictions in real time, training and prediction times are tracked in order to determine the most suitable solution also in terms of responsiveness. In Table 1 it is shown the RMSE at 30 minute and at 60 for both the trained models together with the average time necessary to train the models. As expected, the algorithms are more efficient in the case of shorter time horizons. 5. EXPERIMENTAL RESULTS
Compared to the ARIMA-based model, in general, the LSTM-based model is slightly less performing in terms of RMSE. This, as mentioned, is mainly due to the fact that the ARIMA-based model adapt at every moment in time to the historical series up to that current instant by extracting short or very short-term dependencies, while the LSTMbased model are static models, i.e. they do not change with each new detection, but try to extract longer-term dependencies. Moreover, for the evaluation of the model with a view to industrialization of the solution, as mentioned above, the training time is also relevant. Indeed, for each patient and for each new data point (which is estimated to be collected approximately every 3 minutes using the Glunovo CGM tool) the model is retrained considering the historical and the new detection in order to make a prediction in the next 30 or 60 minutes. For ARIMA, the average training time varies as the parameters p, d and q selected by the developed algorithm vary, but on average the measured value is 0.5 seconds, while for the LSTM the average training time si of 106 seconds. This means that the solution based on LSTM does not allow for real time training and prediction in an online learning fashion. Thus, to use that solution it would be necessary to train a model for each patient periodically, store it and use that model during prediction. This implies that a pipeline to manage models update is necessary. The solution based on ARIMA, on the other hand, allows not to store the model, but just to train it and to discard it after the prediction have been performed. Thus, the LSTM-based solution poses a major scalability problem, since as the number of patients increases, the number of models to re-train and manage will increase, resulting in increased storage and training costs. Therefore, even if the model is good in terms of performance, its industrialization is complex to maintain. The solution based on ARIMA, instead, allows to implement the model in an industrial environment, by training the model with each new collected data for each patient, guaranteeing the prediction in real time.
This study is carried out with the aim of being used by the pilot company within an
application program developed on Google Cloud Platform that continuously monitors
some vital parameters (such as glycemic index, the lipid profile and HbA1c, the blood
pressure and heart rate) and that aims at predicting the trend of the glycemic index in
order to raise alerts. The architecture of the machine learning module implemented is
described in the following. Firstly, the data coming from the sensors are collected and
sent to the time series forecasting module via Google Cloud PubSub [
This research presents a comparison of a statistical and a deep learning based approach for the implementation of a glycemic index time series prediction, monitored by wearable devices. We compare two state of the art approaches: one based on the implementation of statistical ARIMA models by designing an algorithm to automatically identify the value of the optimal parameters of the model based on the use of the ACF and PACF, so that at each iteration the most accurate model possible is created in an online fashion; one based on LSTM with the aim of creating a single prediction model for each monitored patient by taking longer histories into consideration, providing for periodic offline retraining. We used a public dataset to train and test these models which contained real blood glucose monitoring data for different patients. Tests were conducted considering the RMSE with a 30-minute and 60-minute horizon and the time necessary to train and optimize the algorithm. The training was done considering the series in the previous 8 hours. The model based on ARIMA, with our custom parameters optimization algorithm, is the most performing in the 30 and 60 minutes with a RMSE 30 of 1.11 mmol/l and a RMSE 60 of 2.23 mmol/l. Additionally, this model has an average training time of 0.5 seconds which makes it a good candidate for industrialization. On the other hand, if the LSTM-based model has good performances in terms of RMSE (albeit lower than that of the ARIMA-based model), the serving logic of these models is complex: it would be necessary to train a model for each patient with the consequent explosion of the resources necessary for saving, serving and periodic re-training of these models. Thus, we show that the ARIMA-based model can be used in a production environment in an online learning fashion guaranteeing good performances. However, some enhancement can be applied. Indeed, recent studies show how it is possible to increase the performance of glucose level prediction algorithms by considering further data as input, such as information about meals consumed by patients, the doses of insulin taken and the physical exercise performed. All this information can thus be used to improve time series prediction performance. Moreover, in future study, we intend to test the proposed auto adaptive algorithm also in other time series forecasting contexts, like sales forecasting, to test its generalizability. 6. CONCLUSION AND FUTURE WORKS