Improved forecasting of blood glucose could aid in the management of diabetes. Recently proposed neural network architectures that use stacked fully connected layers with residual backcasting have achieved state-of-the-art performance on benchmark timeseries forecasting tasks. Though promising, previous work ignores opportunities for additional supervision, and the use of fully connected layers fails to account for the temporal nature of the signal. Here, we propose a new architecture that builds on previous work by learning to forecast gradually in stages or blocks. Our updates include replacing the fully connected block structure with a recurrent neural network and adding additional losses to provide auxiliary supervision. In addition, we leverage important context in the form of additional input signals. Applied to the task of glucose forecasting, we find that each of these modifications offers an improvement in performance. On the task of predicting blood glucose values 30 minutes into the future, our proposed approach achieves a mean rMSE of 18.2 mg/dL, versus 21.2 mg/dL achieved by the baseline. These improvements get us closer to predictions that are reliable enough to be used in CGMs or insulin pumps to manage diabetes.
Accurate blood glucose forecasting would improve diabetes
treatment by enabling proactive treatment [
We replace the fully connected architecture in each block with a
recurrent neural network (RNN) [
Glucose Additional Var. Forecast
Example Block 1
b1 Block 2 Block 3 Block i
LSTM Fully Connected Split Output bi fi Reconstruction Loss
Forecast Loss
Magnitude Loss
We include additional variables as input to the model (e.g., bolus insulin data), which we hypothesize will provide valuable context. However, due to the sparse nature of some of these variables, we backcast only on our primary variable of interest: blood glucose. We include additional loss terms that act as auxiliary supervision. We hypothesize that these terms will further encourage the model’s blocks to gradually learn accurate components of the signal (i.e., backcasts and forecasts) that sum to the correct signal.
Additionally, we study the effects of pre-training with related datasets. Applied to the task of glucose forecasting, these modifications lead to notable improvements over baseline. 2
Given the strong performance of deep residual forecasting across a
range of tasks [
We focus on the task of univariate time-series forecasting in which we aim to predict future values of a single variable, but assume we have access to additional inputs. Let X = [x(1); :::; x(d)] represent a multivariate time series where d is the number of variables. For each variable i: x(i) 2 RT is a sequence of length T . For our problem, x(1) corresponds to glucose measurements and X0 = [x(2); :::; x(d)] represents other variables of interest (e.g., bolus insulin). Given X, we aim to predict the next h glucose measurements y = x(T1+)1; x(1)
T +2; :::; x(T1+)h. 2.2
We build on a recent univariate time-series forecasting architecture,
N-BEATS [
i=1 fi. The residual nature of the prediction
means that each block learns only what the previous block could
not. Thus, the model learns to gradually reconstruct the signal, while
predicting components of the forecast. We build on this idea by
identifying several opportunities for improvement. Our modification
are as follows:
Accounting for Temporal Structure. We account for the temporal
nature of the input sequence by using a bidirectional LSTM network
[
variables as input at each block. This provides important context when backcasting and forecasting the main signal of interest. However, rather than simply augmenting the dimensionality of the signal throughout, we backcast/forecast only the primary signal, i.e., each block’s input is [xi(1)1 bi 1; X0] (Figure 1). Backcasting/forecasting was not performed on additional variables due to their sparse nature. We hypothesized that learning to backcast/forecast abrupt carbohydrate and bolus inputs would increase the difficulty of the overall task and decrease performance on the main task.
on the loss of the final forecast prediction, M SE(y; y^) = h1 Pih=1(yi y^i)2, we include additional losses with the goal of improving the quality of intermediate representations. We add three auxiliary tasks to our model (shown in Figure 1):
Per-block reconstruction loss (Lr): We hypothesized that explicitly supervising the backcasts could improve their accuracy, and supervising at the block level could improve the intermediate reconstructions. To do this, we add an MSE loss, M SE(xi(1); bi) after each block’s output. Each block’s loss is weighted proportional to its position in the network, encouraging the network to gradually construct a backcast. The final loss Lr is scaled by the sum of the weights of each block, and summed over n blocks: Lr =
i=1 i M SE(xi(1); bi) Pn
Pn i=1 i Per-block forecast loss (Lf ): To encourage each block to contribute towards an accurate forecast, we calculate the loss on the running forecast after each block. We apply a similar per-block scaling term as in the backcast, but increase the weight using cubic values in order to emphasize the network’s later predictions. Let f 0i denote the sum of all forecasts up to the ith block, ie fi0 = Pik=1 fk. The total forecast loss Lf is:
Lf =
Pn i=1 i 3
M SE(fi0; y) Pn i=1 i3 Per-block magnitude loss (Lm): Intuitively, deep residual forecasting aims to incrementally learn parts of a signal. To ensure that each block contributes to the backcast, we penalize blocks by the inverse of their output size (measured using an L1 norm) and scale this penalty by the inverse of their position in the stack, so earlier blocks are encouraged to output larger values. We apply this loss to only the backcast. We do not penalize low-magnitude forecasts because it is possible that a block could account for part of the input signal that is associated with a zero-value forecast. The complete magnitude loss is:
Lm =
Pn 1 1 i=1 i jbij Pn 1 i=1 i Loss is calculated as a weighted sum of these three losses: L = Lf +
Lr +
Lm; where , 2 R+are hyperparameters. Note that the final forecast is included as a part of Lf . 2.3
We ensemble models that use different input lengths. As others have
shown, this allows us to account for trends at different time scales
[
Datasets. We evaluate the proposed approach using the 2020
OhioT1DM dataset [
N-BEATS and uses hyperparameters similar to those reported by
Oreshkin et al.: 10 blocks with 4 layers each, with 512 hidden units
output by each layer, with a ReLU activation function between each
layer [
parameters of our proposed approach, we use Adam [
To reduce the chance of overfitting to the 2020 data, we tune all hyperparameters, excluding and , on the six subjects from the previous competition. We performed a grid search to select the number of blocks (range: 5-10, best: 7) and number of hidden LSTM units (range: 50-300, best: 300), optimizing for rMSE of the final time step of the prediction on the held out data. Including more hidden units could perhaps have further improved performance, but we were limited by memory constraints.
During tuning, we measured the added value of including each auxiliary variable (e.g., sleep data), including only those variables that improved performance for the final evaluation. We tuned the block-weighting terms for our auxiliary loss functions, considering inverse, constant, and linear relationships for Lr, and powers of 2 and 3 for Lf . We tuned the auxiliary loss function weights , and using the validation data for the new subjects. We performed a grid search over [:1; :15; :25; :3; :5; 1; 2; 4] for and 10e5 times those values for , selecting = :3 and = 10e4. Grid search analyses were tuned with a subset of backcast lengths (3f and 6f ) across all patients, for efficiency.
For both the 30 and 60 minute prediction horizons, we train using loss for the full window (i.e., h = 6 for the 30 minute analyses and h = 12 for the hour prediction), but report results using the last time step only (as dictated by the competition).
Pre-training. As mentioned above, we pre-train on the 2018
OhioT1DM and Tidepool data. Pooling these datasets together,
we train a single model. The weights of this model are used to
initialize the weights of a ‘global’ model trained on the 2020
OhioT1DM training set. This ‘global’ model is then fine-tuned to
each participant from the 2020 OhioT1DM dataset, resulting in six
participant-specific models. Baseline analyses do not use the large
pooled dataset for pre-training, but start from a single global model
(trained from random initialization) and then fine-tune.
Code. All of our experiments were implemented in Python/PyTorch
[
We report results for models trained to predict 12 time steps (one hour horizon) and six time steps (30 minute horizon), reporting outcomes on only the final time step. We evaluate using the square root of the mean squared error (rMSE), and mean absolute error (MAE), defined as M AE(y; y’) = jy ny0j for n samples. Metrics are applied directly to the raw CGM values with no preprocessing.
To examine the effect of each architecture modification, we perform an ablation study where we add each modification to the architecture incrementally (i.e., first RNN blocks only, then RNN blocks with added variables, then RNN blocks with added variables and additional loss terms, then the full architecture with pre-training). We also experiment with each modification made in isolation.
Finally, beyond the evaluation dictated by the competition, we
also calculate outcomes for the most crucial forecasting windows:
those that represent the beginning of a hyper- or hypoglycemic event.
These correspond to predictions for which the most recent time step
is in the euglycemic range (blood glucose from 70-180 mg/dL),
but the participant becomes hypoglycemic (< 70 mg/dL) or
hyperglycemic (> 180 mg/dL) within the prediction horizon. We also
examine the distribution of predictions vs. actual measurements in
terms of the corresponding would-be treatment recommendations
with a Clarke Error Grid Analysis [
Overall, the proposed approach leads to average rMSEs of 18.2 and
31.7 on the 30 and 60 minute horizons, respectively, and average
2 https://gitlab.eecs.umich.edu/mld3/deep-residual-time-series-forecasting
MAEs of 12.8 and 23.6 for those horizons (Table 1). These values
are comparable to the results achieved in the 2018 challenge (i.e.,
18.9 to 21.7 for 30 minute rMSE) [
Based on our ablation study, each of our modifications improves performance over baseline, although to varying degrees (Figure 2). For the 30 minute horizon, mean rMSE across participants is 21.2 for the baseline model (BL). Adding the RNN block structure (bidirectional LSTM) and additional variables (+AV) improves performance to 20.1 and then 18.94, respectively. Adding the additional loss functions (+LOSSES) results in a slightly improved mean performance of 18.87, and pre-training on the Tidepool and 2018 OhioT1DM datasets (+PT) further improves performance to 18.2. We observed similar trends for the 60 minute horizon, and for MAE.
While all of our modifications reduce rMSE, one of the largest improvements in performance comes from replacing the fully connected layers of N-BEATS with recurrent layers (improves performance from 21.2 to 20.1). These recurrent layers directly model each time step as a function of the previous, allowing for more accurate temporal representations. Compared to the baseline, the proposed approach leads to more faithful predictions and fewer extreme failures (Figure 3). Even without pre-training, the modified architecture eliminates the extreme failures observed in the baseline case (highlighted in Figure 3).
Additional variables are significantly more impactful when added to the RNN-based block structure (a 30 minute rMSE improvement from 20.1 to 18.9), when compared to adding them in isolation (an improvement from 21.2 to 21.0). Despite the poor performance when used in the original baseline architecture, we chose to include these variables in our updated RNN-based model. The RNN block structure has multiple input channels, which allows for the direct use of the variables. In contrast, the fully connected layers require a flat input, which removes much of the timing information. This illustrates the importance of the ordering of decision making in model architecture development. Although such decisions are easy when an exhaustive search of combinations is performed, this is not always feasible, so these decisions often must be made sequentially. If we had selected additional variables based on the performance of the baseline model alone, we may not have opted to use them, and would have lost the large benefit that they provide when used in tandem with the RNN block structure.
To test our hypothesis regarding the efficacy of backcasting on blood glucose values only (and not additional variables), we ran a post-hoc analysis with a model using residual backcasting for all input variables. As expected, we found that this decreased performance (30 minute rMSE = 19.2 vs. 18.9).
In contrast to the other modifications, the additional loss terms offer only slight improvements over the baseline. During model selection, this modification offered substantially more improvement on the OhioT1DM 2018 dataset (decreasing rMSE of the RNN + added variables model from 19.6 to 18.9). Yet even with additional tuning, this modification is only minimally helpful on the 2020 dataset, suggesting that the method may not be universally beneficial. In an ablation analysis, a model using only (Lf ) performs best (rMSE = 18.82) on the 2020 data. However, in contrast, the added variables seem to help more for the 2020 data than for the 2018 data. This suggests that perhaps the losses serve a regularizing effect when insufficient information is present in the additional variables, as we observed in the 2018 validation data.
Pre-training improves the performance of both the modified and baseline architectures, although it offers the most improvement for the baseline (improving performance from 21.2 to 19.8 when applied alone). This suggests that pre-training becomes less crucial once other improvements are made. Running an analysis using our proposed architecture and no pretraining over 4 random seeds results in similar performance across seeds (30-minute rMSEs were 18.87, 18.82, 18.95, and 18.89), indicating fairly robust performance over different random initializations. Further, pre-training improves performance significantly over all random initialization, demonstrating the benefit of pre-training beyond chance.
In our event-specific evaluation (includes only prediction windows where a hyperglycemic or hypoglycemic event occur), the mean rMSE is significantly higher compared to the holistic evaluation (for our proposed approach 30 minute horizon: 27.9 vs. 18.2). This is expected, since such events are difficult to anticipate. Specifically, the proposed approach achieves a higher rMSE for hyperglycemic events (30.5) than hypoglycemic (23.8). Interestingly, while our model performs better than the baseline for all events on average (our model: 27.9 vs baseline: 31.0) and for hyperglycemic events (30.5 vs 34.2), the two models perform comparably in the hypoglycemic range (23.8 vs 23.6). This could be due to the relative rarity of hypoglycemic events (for example, there were 322 in the test data, vs 868 hyperbe recommended based on our proposed approach’s predictions vs what treatments should be given based on the actual value. 99% of predictions fall within regions A and B, which indicate appropriate treatment recommendations. (b) Proportion of predictions in each region of the grid for our proposed approach and the baseline model. glycemic events), or it could be a reflection of the MSE loss function over-emphasizing large values, which is unaccounted for in our approach.
In our Clarke Error Grid Analysis [
We find deep residual forecasting to be effective when applied to the task of predicting blood glucose values. Augmenting a previously proposed architecture with RNNs in place of fully connected stacks, additional variables, and self-supervising loss functions all lead to improvements when applied to the task of blood glucose forecasting. Beyond blood glucose forecasting, we hypothesize that many of the proposed changes could be beneficial when applied to other forecasting tasks.
This work was supported by JDRF (award no. 1-SRA-2019-824-SB). In particular, this award provided access to the Tidepool dataset that was used in pre-training the models. The views and conclusions in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied of JDRF.