{'fraction_non_alphanumeric': 0.03447446002045605, 'fraction_numerical': 0.009626376270982492, 'mean_word_length': 4.8985095812633075, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gluformer: Transformer-Based Personalized Glucose Forecasting with Uncertainty Quantification

Deep learning models achieve state-of-the art results in predicting blood glucose trajectories, with a wide range of architectures being proposed. However, the adaptation of such models in clinical practice is slow, largely due to the lack of uncertainty quantification of provided predictions. In this work, we propose to model the future glucose trajectory conditioned on the past as an infinite mixture of basis distributions (i.e., Gaussian, Laplace, etc.). This change allows us to learn the uncertainty and predict more accurately in the cases when the trajectory has a heterogeneous or multi-modal distribution. To estimate the parameters of the predictive distribution, we utilize the Transformer architecture. We empirically demonstrate the superiority of our method over existing state-of-the-art techniques both in terms of accuracy and uncertainty on the synthetic and benchmark glucose data sets.


INTRODUCTION
Prediction of blood glucose values in patients with diabetes is an active area of research [1]. From a clinical standpoint, accurate forecasting of glucose levels can help patients take proactive actions and prevent severe complications such as hypoglycemia or diabetic coma [2,3]. From the methodological perspective, the highly non-linear and non-stationary nature of the glucose profiles makes accurate predictions difficult. To the best of our knowledge, most current methods produce point-wise predictions without quantifying their uncertainty, which limits clinical applicability. A predicted glucose increase may prompt the patient to take extra insulin, which may lead to dangerous hypoglycemia if the prediction is wrong. Complementing predictions with uncertainty quantification will help to minimize these risks.
Since existing approaches primarily focus on maximizing the accuracy of the predictions (rather than quantifying uncertainty), they focus on minimizing the mean squared error (MSE) between the predictions and the observed data [4,5,6,7,8,9,10,11,12]. MSE loss implicitly corresponds to modeling future trajectory as a function of the past plus a homogeneous Gaussian noise. If the assumed model is correct, this approach leads to the best predictions, and confidence intervals can be constructed by inferring the corresponding noise variance. However, a homogeneous Gaussian noise assumption is typically violated in practice, making the corresponding predictions and confidence intervals unreliable.
Here we propose a new framework for modeling future glucose trajectories using an infinite mixture model (IMM). Such formulation allows us to capture a variety of non-Gaussian output distributions, e.g., multi-modal, skewed. We estimate the parameters of the IMM using the Transformer architecture. While IMM has an infinite number of parameters, we introduce an approximation step to fit the distribution implicitly. Specifically, we inject randomness into the computation of the network outputs and treat each output as a realization from one of the mixture components.
To validate the capabilities of the proposed model, we consider two data sets: 1) synthetic data set 2) publicly available continuous glucose monitor (CGM) data set [11]. We achieve state-of-the-art results for both. Figure 1 shows a sample of predicted curves together with the visualization of predictive distribution, demonstrating the superiority of our method in terms of the robustness of uncertainty estimates to deviations from the Gaussian assumption.

RELATED WORK
Traditional machine learning. ARIMA is a well-known model for time series data and is commonly used with CGM data [4,5,7,8]. ARIMA is a non-linear auto-regressive approach that extends the standard AR model to handle nonstationarity by integrating the previous trends. Another popular model for glucose predictions is random forest [6,9,10]. The model splits the data into random subsamples, fitting a single prediction tree for each. The average of the ensemble of individual trees is used for prediction. The model is typically applied recursively by producing a single prediction at a time, but multi-output extensions are also available. Deep learning. As time-series data closely resemble data from the natural language processing field, there has been an interest in applying deep learning for multi-output timeseries forecasting. Some existing state-of-the-art approaches for glucose prediction are PolySeqMO [11] and RNN [12]. PolySeqMo models the predicted curve in terms of a polynomial expansion and uses an RNN to learn the coefficients of the polynomial. The RNN approach of [12] improves Poly-SeqMO by moving away from the polynomial expansion and building upon the approach taken in [13], providing personalized predictions, proposing attention for longer inputs, and introducing robust training scheme for the model. Recently, the Temporal Fusion Transformer (TFT) has shown promising results on a variety of time series data sets [14]. The model is formulated using quantile regression and hence allows for uncertainty quantification. However, unlike our model, the TFT does not provide any measure of the distributional fit, which prevents the users from judging the quality of the TFT uncertainty estimates. Uncertainty estimation. The infinite mixture model (IMM) has been previously studied in [15,16,17]. These approaches rely on MCMC sampling from the posterior distribution, with latent variables explicitly being sampled from specified prior. In contrast, we estimate the IMM parameters implicitly by injecting latent variables (noise) into the network. This change allows us to translate the sampling problem into the maximum-likelihood optimization, where the only additional cost we pay is several stochastic forward passes through the network. Since full MCMC sampling is costly to implement and diagnose for high-dimensional data sets, our approach yields a significant computational advantage.
Stochastic neural networks have been investigate by [18] and [19]. [18] explore a connection between a stochastic neural network with a dropout and L 2 regularization and the Gaussian process. Specifically, they show that the dropout imposes a distribution on the weights of the neural network, which can then be optimized as a variational approximation to the Gaussian kernel parameters. Their derivations can be seen as a special case of our model if we take the basis distribution to be Gaussian, impose L 2 constraint on the weights, and implement the latent variables as a dropout.

Notation and background
The blood glucose prediction can be viewed as a time series forecasting problem. Let x i = (x ij ) t j=1 be a sequence of glucose values of length t with i = {1, . . . , n} being the sample index. Denote the consecutive future glucose values as . . , y iT ) forms an equally-spaced time sequence of CGM readings. Each such sequence corresponds to a unique patient in the data set. The goal is to estimate y given history x.
To achieve this goal, previous works (Section 2) focus on building a predictive model f θ : R t → R T with parameters θ. These models usually rely on the parametric assumption y ∼ N (f θ (x), σ 2 I), which is equivalent to the additive Gaussian assumption on the residuals. While this assumption is sometimes explicit [13], it is often implicit through the choice of the loss function. When the true predictive distribution, p(y|x), is far from the Gaussian, such models are clearly misspecified. Subsequently, both predictions and uncertainty estimates of such models become unreliable.

The Infinite Mixture Model
We first consider univariate prediction and denote the scalar response variable as y. We introduce latent (hidden) variables z ∼ q, where q is some distribution that allows fast sampling. Provided we have sampled z from q, we assume that the base distribution assumption holds conditionally on z: where s is a vector of the sufficient statistics for a distribution p, e.g., mean and variance for Gaussian. Marginalizing over z, we get the full predictive distribution expressed as an infinite mixture model (IMM): p(y|x) = p(y; s(x, z))q(z)dz We propose to estimate the sufficient statistics s with a neural network modelŝ := f θ (x|z). To approximate the intractable predictive distribution, we use Monte Carlo integration. That is, we sample k latent z j independently from q, and average over k stochastic passes through the network: The latent distribution, q, can be implemented by adding a random noise or concatenating a latent noise vector to the input vector, x. Alternatively, the latent distribution may be embedded in the network structure through the Bernoulli dropout layers. The choice of the base distribution, p, ultimately depends on the specifics of the problem. While the Gaussian base distribution is a good starting point due to its theoretical properties [15], Laplace base distribution may be preferable when modeling highly asymmetric distributions [20].
To generalize from the univariate case to the multi-output forecasting problem, we model the multivariate predictive distribution p(y|x) by specifying univariate IMM models as above for each output and combining them using an appropriate copula dependence function (e.g., independence copula, Gaussian copula etc.). Since, by Sklar's theorem, any multivariate distribution can be decomposed into its marginals and a dependence function, our specification remains expressive [21]. The copula function parameters can then also be estimated using the neural network. In practice, we found that the independence copula, corresponding to the independence across time assumption, was sufficient.

Measuring the Quality of Predictive Distribution Fit
Likelihood. To assess the quality of estimated predictive distributions, we compute log-likelihoods for each method on the test data, where a higher value indicates a better fit.
For existing models (Section 2), the additive Gaussian assumption on the noise leads to the predictive distribution y|x ∼ N (f θ (x), σ 2 I). Since such models do not explicitly estimate the variance, σ 2 , during fitting, we estimate it using MLE and compute the average log-likelihood as: For our model, we approximate the average log-likelihood on the test data set based on the MC draws. That is, we estimate the average log-likelihood as: Calibration. To further assess the quality of our predictive distribution on the benchmark CGM data set, we propose to evaluate model calibration. We focus on the quality of the estimated marginal distributions for each predicted time point. Given n input sequences, our approach produces approximate predictive distributions based on the MC draws, {p(y i |x i )} n i=1 . Based on our copula formulation for multivariate output, we also have direct access to the marginal distributions, {[p(y ij |x i )] T j=1 } n i=1 , each of which is modeled as an IMM. Thus, we can derive approximate marginal cumulative density functions (CDF) as: where F is the CDF of the chosen univariate base distribution.
To judge the calibration of our model (how close are the estimated confidence intervals to nominal coverage) via a diagnostic plot, we utilize the approach described in [22]. We select a set of η jl ∈ [0, 1] for each j ∈ {1, . . . , T } and a pre-specified mesh l ∈ {1, . . . , L}, and calculatê Thus, we obtain a set of points {(η jl , η jl )} j for each l ∈ {1, . . . , L}, whereη jl are empirical frequencies and η jl are expected (true) frequencies. Ideally, in the case of perfect calibration, we have thatη jl = η jl .

Estimation
Transformer Architecture. The Transformer model [23] consists of a series of attention heads and fully connected layers. To overcome the computational cost of modeling longer prediction intervals, we use a modification of the Transformer architecture that introduces convolutional layers in-between attention blocks and is able to provide multioutput predictions in a single pass [24]. We use the model to predict the parameters of each of the marginal IMM's, where we consider each stochastic pass through the model to be an MC draw. Feature Embedding. As outlined in Section 4.1, the CGM data usually only contains glucose readings together with the timestamps. It is reasonable to assume that each subject adheres to a daily schedule that affects his or her blood glucose levels. Therefore, we utilize time features to enhance our predictions. Also, we add positional and subject encoding.  Table 2. APE/RMSE for 60-minute prediction window (Full), hypoglycemia, hyperglycemia, and event (hypo-or hyperglycemia), and model log-likelihood on test data.  Figure 1. Benchmark CGM. We use a publicly available CGM data set [11], which contains information on glucose levels of 38 subjects tracked continuously throughout multiple disjoint intervals with the measurement frequency of 5 minutes. Similar to [12], we remove periods of drastic fluctuations where the subsequent CGM readings change by more than 40 mg/dL. To be compatible with the analysis in [12], we do not utilize any interpolation techniques for the missing data. The resulting data consists of 399, 302 observations, with an average of 30 (uninterrupted) sequences per subject of length between 200 and 400 observations. We split the data by sequences into the train, validation, and test sets in 20 : 1 : 1 proportion.

Model Implementation
We implement our method in PyTorch [25], and use a single NVIDIA RTX 2080 Ti GPU for all experiments. For ablation purposes, we train the model using (i) proposed loss corresponding to the IMM; (ii) mean squared error (MSE) corresponding to the Gaussian assumption. For IMM, we implement the latent variables by introducing the dropout and assuming independence across time. The full model implementation is available at github.com/gluformer. For the synthetic data, we set the encoder length to 4 and the prediction to 2. For the glucose data, we set the encoder length to 180 (36 hours) and the prediction length to 12 (1 hour).

Results
We measure prediction accuracy using root mean squared error (RMSE) and average percentage error (APE) across prediction intervals. We report the final APE and RMSE as medians of the individual prediction errors across samples. We measure the quality of fit as in Section 3.3. Synthetic. We compare the proposed IMM to the same architecture trained with the Gaussian assumption. We report our results in Table 1 and plot a sample of predictions in Figure 1, displaying the predicted densities for the two models. The proposed IMM is more accurate and correctly captures the multi-modality in the predicted densities across time points, providing a better probabilistic fit to the data. Benchmark CGM. We report both the error on the full data set as well as errors for the periods of hypo-and hyperglycemia [26]. We define a period of hypoglycemia if the glucose level drops below 70 mg/dL in the forecast window and a period of hyperglycemia if the glucose level rises above 180 mg/dL. We also measure the model's accuracy during the incidence of either of the events. Table 2 summarize the results on the benchmark CGM data for the 60 minutes prediction window. The proposed IMM outperforms the others 3 out of 4 times in APE and 2 out of 4 times in RMSE. In particular, IMM achieves state-of-the-art 7.78 APE for the most difficult 60 minute forecast. Table 2 also reports average likelihoods on test data for all methods, except TFT for which it is not available. Our model achieves 6× better average log-likelihood than the previous methods. We believe this to be indicative of the fact that the true predictive distribution of y|x is not similar to Gaussian (likely multi-modal). We report our calibration plots for each of the 12 prediction points in the future (the longest prediction length) in Figure 2. When plotted, a perfectly calibrated model should correspond to a straight line. We find that the model observed confidence levels are close to the expected levels for all the points, highlighting the model's excellent performance.

CONCLUSION
We propose a novel approach to quantify uncertainty in multioutput time series predictions based on the infinite mixture model. We utilize an efficient Transformer architecture for estimation. When applied to public CGM data, our model outperforms existing approaches (ARIMA random forest, RNN, and Transformer) in terms of the achieved log-likelihood on the test data and provides competitive accuracy for various forecast horizons and dangerous event scenarios.