Kernel machine learning methods to handle missing responses with complex predictors: application in modelling five-year glucose changes using distributional representations

Computer Methods and Programs in Biomedicine 221 (2022) 106905 Contents lists available at ScienceDirect Computer Methods and Programs in Biomedicine journal homepage: www.elsevier.com/locate/cmpb Kernel machine learning methods to handle missing responses with complex predictors. Application in modelling five-year glucose changes using distributional representations Marcos Matabuena a , ∗, Paulo Félix a , Carlos García-Meixide b , Francisco Gude c a CiTIUS (Centro Singular de Investigación en Tecnoloxías Intelixentes), Universidade de Santiago de Compostela, Santiago de Compostela 15782, Spain b ETH Zürich c Unidade de Epidemioloxía Clínica, Complexo Hospitalario Universidade de Santiago (CHUS), Travesía da Choupana, Santiago de Compostela 15706, Spain a r t i c l e i n f o Article history: Received 14 October 2021 Revised 11 May 2022 Accepted 22 May 2022 Keywords: Missing data Kernel methods Statistical independence Variable selection Regression modelling Diabetes mellitus Continuous glucose monitoring a b s t r a c t Background and objectives: Missing data is a ubiquitous problem in longitudinal studies due to the number of patients lost to follow-up. Kernel methods have enriched the machine learning field by successfully managing non-vectorial predictors, such as graphs, strings, and probability distributions, and have emerged as a promising tool for the analysis of complex data stemming from modern healthcare. This paper proposes a new set of kernel methods to handle missing data in the response variables. These methods will be applied to predict long-term changes in glycated haemoglobin (A1c), the primary biomarker used to diagnose and monitor the progression of diabetes mellitus, making emphasis on exploring the predictive potential of continuous glucose monitoring (CGM). Methods: We propose a new framework of non-linear kernel methods for testing statistical independence, selecting relevant predictors, and quantifying the uncertainty of the resultant predictive models. As a novelty in the clinical analysis, we used a distributional representation of CGM as a predictor and compared its performance with that of traditional diabetes biomarkers. Results: The results show that, after the incorporation of CGM information, predictive ability increases from R 2 = 0 . 61 to R 2 = 0 . 71 . In addition, uncertainty analysis is useful for characterising some subpopulations where predictivity is worsened, and a more personalised clinical follow-up is advisable according to expected patient uncertainty in glucose values. Conclusions: The proposed methods have proven to deal effectively with missing data. They also have the potential to improve the results of predictive tasks by including new complex objects as explanatory variables and modelling arbitrary dependence relations. The application of these methods to a longitudinal study of diabetes showed that the inclusion of a distributional representation of CGM data provides greater sensitivity in predicting five-year A1c changes than classical diabetes biomarkers and traditional CGM metrics. ©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) 1. Introduction Missing data are common in epidemiological studies. On the face of it, the extended practice of excluding participants with only partially available data on the variables of interest results in ignoring valuable information, thereby leading to biased estimates which often rely on unrealistic assumptions [1–3] . To draw reli- ∗Corresponding author. E-mail address: mar[email protected] (M. Matabuena) . able conclusions, principled methods are imperative by appropriate modeling of the missingness mechanism [4] . Kernel methods are a class of effective pattern recognition algorithms that are well suited to model nonlinear relations between the response and predictors. These are built on the notion of a kernel function as a similarity function between a new instance and those included in the training set [5,6] . One of the most significant achievements of kernel methods is the proposal of appropriate kernel functions for managing complex statistical objects such as graphs, strings, or probability distributions [7] . Thus, kernel methods are expanding the range of possible applications for https://doi.org/10.1016/j.cmpb.2022.106905 0169-2607/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) M. Matabuena, P. Félix, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 221 (2022) 106905 Fig. 1. Left: The 5-day CGM recording from a normoglycemic patient is shown. Centre: Glucodensity designates a distributional representation that estimates the proportion of time the patient spent at each glucose concentration. Right: Quantile representation. Dotted, solid and dashed lines represent concentrations for 20 percent, 50 percent and 80 percent quantiles, respectively. machine learning in the health domain, challenged with the rapid increase in new complex medical data. The main purpose of this study is to propose a set of kernel methods ( Section 2 ) to handle missing responses for statistical independence testing ( Section 2.2 ), variable selection ( Section 2.3 ), and conformal inference ( Section 2.4 ). One major advantage of these methods is their ability to operate as a sequence of predictive stages which increasingly filter out irrelevant information, while also providing an evaluation of the limits of the ensuing predictions. In particular, the present proposal is based on the reproducing kernel Hilbert space (RKHS) 1 framework, providing a Hilbert space of functions that is fully characterized by a reproducing kernel. Importantly, every function in an RKHS that minimizes an empirical risk function can be written as a linear combination of the kernel function evaluated at the training data, and it is ensured that a solution for a machine learning problem that is close to the true solution and also generalizes well to the test data can be obtained. An essential property of the RKHS framework is that it overcomes the limitations of previous proposals focused on Euclidean functional representations [8] . The proposed methods are motivated by the need to explore the limits of predicting long-term glucose changes in a five-year longitudinal population-based study, including both healthy and diabetic individuals, where a subsample of participants underwent continuous glucose monitoring procedures at the beginning of the study. As expected, a substantial number of participants withdrew from the study, and therefore, an analysis robust to missing values in the response variable is required in order to maintain the validity of the statistical inferences [9] . We include a novel distributional representation for CGM data as a predictor (see Fig. 1 ) [10] . Among the different biomarkers, we select the glycated haemoglobin (A1c) as the response variable. A1c is a measure of the average blood glucose level over the past three months, and it is the preferred option because it provides more reproducible values in the laboratory and is subject to less measurement error [11] . Furthermore, we aimed to assess and discuss the residuals and predictive capacity of several variables associated with the evolution of A1c in the long term, providing interpretable clinical phenotypes for large uncertainty cases. 1 Appendix contains a guide to acronyms used in this document. The rest of this paper is outlined as follows. Section 2 describes in detail the methods for statistical independence testing, variable selection, and inference on the uncertainty of new predictions. Section 3 presents the application of these methods for modeling long-term changes in glucose. Section 3.1 describes the AEGIS database used to test the proposed method. Section 4 presents the results of the application of these methods to the AEGIS database. Section 5 discusses the advantages and drawbacks of the proposed approach. Finally, conclusions are presented in Section 6 . 1.1. Data analysis outline Finally, this paper presents a data analysis framework designed as a pipeline of kernel methods for predictive problems with missing data. Their subsequent application to diabetes mellitus will allow us to examine the relationship between the baseline characteristics of participants in a five-year study and A1c as the response variable. The proposed framework comprises the following steps: 1. To measure the statistical association between each predictor and the response variable with an efficient statistical independence test. If the response is proven to be independent of a predictor, it can be screened out from further consideration. To this end, we adapted a previous kernel independence test and designed a new bootstrap method to perform test calibration. The test was applied to check the association between some diabetes biomarkers and five-year changes in the A1c variable, A1c 5 years −A1c initial . 2. To identify the best subset of predictors revealing higher-order interactions with the response variable in order to improve the prediction. To this end, we adapted a previous kernel variable selection method and applied it to find the best subset of diabetes biomarkers most strongly associated with A1c 5 years . 3. To explore the prediction ability of a set of explanatory variables through a non-linear regression method. To this end, we adapted a previous kernel ridge regression method and applied it to predict A1c 5 years . 4. To estimate the uncertainty of the predictions. To this end, we designed a new method to provide a prediction interval for the response variable, based on conformal inference. Using this method, we can measure the limits of the regression models previously obtained and, significantly, identify specific patient 2 M. Matabuena, P. Félix, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 221 (2022) 106905 subpopulations that do not fit the expected behaviour, which is a key issue for clinical decision-making. 2. Methods 2.1. Preliminaries We first pose the problem in general terms. Let (X , Y, R ) ∈ X × R ×{ 0 , 1 } be a random vector such that X = (X 1 , . . . , X p ) denotes the predictor variables, Y is the response variable, and R is a binary random variable that indicates whether the response is missing. X denotes a general topological space, meaning that it can be arbitrary, discrete, continuous, or structured. Let D = { (x i , y i , r i ) } n i =1 be a dataset of independent, identically distributed observations, where y i is missing if r i = 0 . We assume R is distributed according to R | X ∼Ber(π(X)) , with π(·) = P (R = 1 | X = ·) ; hence, some of the predictors can have an impact on the mechanism of missing data π(·) = P (R = 1 | X = ·) . For instance, in our example, older patients are more reluctant to perform a second CGM monitoring, so the probability of not observing a patient increases with age. We also assume a missing at random (MAR) mechanism in the response Y , namely, R and Y are conditionally independent given X or, in short, R ⊥ Y | X . Consider the following relation between X and Y : Y = f(X ) + , (1) where denotes a random noise with E (| X ) = 0 , and fis the true regression function. Our goal is to predict Y by proposing a new data analysis framework that is robust to datasets in which some values for Y are not observed. To this end, we provide: 1) a method for univariate analysis based on testing the statistical independence between each predictor variable and the response variable; 2) a method for selecting the subset of predictor variables that best predicts the response variable; and 3) methods for predicting the response variable and inferring the uncertainty in the predictions. These methods are based on the reproducing kernel Hilbert space (RKHS) learning paradigm. The core element of this paradigm is a positive definite kernel function k X : X ×X → R which allows us to measure the similarity between any x and y , with x, y ∈ X . The positive definiteness of the kernel function guarantees the existence of a dot product space Hand feature mapping φ: X → H, such that k X (x, y ) =  φ(x ) , φ(y )  H . Thus, we can express a broad spectrum of statistical modeling problems as linear [12] , allowing computational algorithms to easily determine optimal solutions. 2.2. Testing statistical independence We wish to test whether two random variables X ∼P X and Y ∼P Y are not independent, that is, if we can reject the null hypotheses H 0 : X ⊥⊥ Y , from n samples { (x i , y i ) } n i =1 . To do this, we must calibrate the test under the null hypothesis to determine the results that are expected to occur with a certain probability if the null hypothesis holds. In our specific case, we must consider the effects of the mechanism of missing data in the response variable Y . We propose a methodology to address this problem based on kernel mean embeddings, which is valid when both covariate and response variables live in a separable Hilbert space. In addition, we introduce a new bootstrap procedure to perform the test calibration adapted to kernel mean embeddings. Hilbert space embeddings of distributions or, in short, kernel mean embeddings [7] , allow us to map distributions into a RKHS, in which kernel methods can be extended to probability measures. Kernel mean embeddings can be used to define a metric for distributions, the maximum mean discrepancy (MMD), which in turn can be applied to define an independence test, the Hilbert-Schmidt independence criterion (HSIC), a non-parametric test of independence with the important property that it does not make any assumption as to the nature of the possible dependence among the two variables [13] . We extended this test to a missing data setting. A reproducing kernel of His a kernel function that satisfies (1) ∀ x ∈ X , k X (·, x ) ∈ H, and (2) ∀ x ∈ X , ∀ f ∈ H,  f, k X (·, x )  H = f(x ) . If Hhas a reproducing kernel, it is said to be an RKHS, H k X . Kernel mean embedding results from extending the mapping φto the space of probability distributions by representing each distribution as a mean function φ(P ) =  X k (·, x ) dP (x ) , resulting in the transformation of a distribution P into an element of the RKHS H k X . Given two probability measures P and Q, the RKHS distance between their embeddings can be defined as the MMD [14] : MMD k X (P, Q) =  φ(P ) −φ(Q)  H k X . (2) For the class of characteristic kernels, the embeddings are injective, that is, MMD k (P, Q) = 0 , if and only if P = Q. MMD can then be applied to measure the degree of dependence between the random variables X ∈ X and Y ∈ Ywith marginal distributions P X and P Y and jointly distributed as P X,Y . Note that testing the null hypothesis H 0 : X ⊥⊥ Y is equivalent to testing H 0 : P X,Y = P X P Y . We denote by φX (·) , φY (·) and φX,Y (·) the kernel mean embeddings of P X , P Y , and P X,Y , respectively. Assuming H k Z is a RKHS over X ×Ywith kernel k Z ((x, y ) , (x  , y  )) = k X (x, x  ) k Y (y, y  ) , so that H k Z is a direct product H k X H k Y (with being the tensor product), then a natural way of testing independence is measuring the MMD distance between the functions φX,Y (·) and φY (·) φX (·) , which can be written as the Hilbert-Schmidt Independence Criterion (HSIC) between Xand Y [14] , defined as HSIC (P X,Y , P X P Y ) = || φX,Y −φX φY || 2 H k X H k Y (3) It can be shown that when k X and k Y are characteristic kernels, HSIC (P X,Y , P X P Y ) = 0 if and only if X ⊥⊥ Y . Expanding Eq. (3) , we have HSIC (P X,Y , P X P Y ) =  φX,Y −φX φY , φX,Y −φX φY  H k X H k Y =  φX,Y , φX,Y  +  φX φY , φX φY  −2  φX,Y , φX φY  , (4) where H k X H k Y is dropped in the subscript for brevity. From the reproducing property, E P [ f(x )] =  f, φ(P )  H , ∀ f ∈ H, and Fubini’s theorem, we obtain HSIC (P X,Y , P X P Y ) = E X,Y,X  ,Y  [ k X (X, X  ) k Y (Y, Y  )] + E X,X  [ k X (X, X  )] E Y,Y  [ k Y (Y, Y  )] − −2 E X,Y [ E X  (k X (X, X  )] E Y  [ k Y (Y, Y  )]] , (5) where X  and Y  are independent copies of random variables X and Y , respectively. Ultimately, testing independence involves calculating the squared distance between two mean functions in the appropriate RKHS space, resulting from transforming the original data to capture all distributional differences between both random variables. In practice, a limited number of samples, { (x i , y i , r i ) } n i =1 , are observed. Therefore, we must replace the population mean with the sample mean, defined through its empirical distribution. Then, the Hilbert-Schmidt independence criterion can be estimated as  HSIC ( ˆ P X,Y , ˆ P X ˆ P Y ) = 1 n 2 n  i =1 n  j=1 k X (x i , x j ) k Y (y i , y j ) + 1 n 2 n  i =1 n  j=1 k X (x i , x j ) n  i =1 1 n 2 n  i =1 n  j=1 k Y (y i , y j ) −1 n 3 n  i =1 n  j=1 n  k =1 k X (x i , x j ) k Y (y i , y k ) . (6) 3 M. Matabuena, P. Félix, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 221 (2022) 106905 Under the MAR assumption, we observe { (x i , y i , r i ) } n i =1 , and we must estimate the missing data mechanism given by the function π(·) = P (R = 1 | X = ·) . Several procedures have been proposed in the literature for this purpose, such as logistic regression, lasso, random forest, and ensemble methods. Subsequently, we re-weighted the dataset, taking into account the difficulty of observing the response of the i th datum. In particular, we associate weight w i with the i th datum via an inverse probability weighting (IPW) estimator [4] given by w i = r i nπ(x i ) , i = 1 , . . . , n, (7) which results in assigning large w i values as the probability of observing a response decreases. Using this procedure, we obtain an asymptotic unbiased estimator that balances the sampling mechanism and allows us to make a proper inference according to the target population examined. We define the normalized weight of w i as w ∗ i = w i  n i =1 w i , i = 1 , . . . , n. (8) We denote the estimated and normalized i th weight as ˆ w i and ˆ w ∗ i , respectively, after estimating ˆ π(·) . To obtain an estimator of HSIC with missing data, it is sufficient to replace the uniform weight 1 /n of the empirical distribution with normalised weights ˆ W ∗= ( ˆ w ∗ 1 , . . . , ˆ w ∗ n ) in Eq. (6) . Thus, we obtain  HSIC ( ˆ P X,Y , ˆ P X ˆ P Y ) = n  i =1 n  j=1 ˆ w ∗ i ˆ w ∗ j k X (x i , x j ) k Y (y i , y j ) + n  i =1 n  j=1 ˆ w ∗ i ˆ w ∗ j k X (x i , x j ) n  i =1 n  j=1 ˆ w ∗ i ˆ w ∗ j k Y (y i , y j ) − n  i =1 n  j=1 n  k =1 ˆ w ∗ i ˆ w ∗ j ˆ w ∗ k k X (x i , x j ) k Y (y i , y k ) . (9) Calibration under the null hypothesis with the precedent statistic is not trivial, and the permutation approach is generally not valid because the response Y is not exchangeable due to the nonhomogeneous missing data mechanism. To overcome this difficulty, we propose a novel bootstrap approach that properly deals with non-vectorial predictors [15] . Under null hypothesis H 0 : P X,Y = P X P Y , it can be assumed that φX,Y (·) −φX (·) φY (·) = 0(·) . Therefore,  HSIC ( ˆ P X,Y , ˆ P X ˆ P Y ) =  ˆ φX,Y −ˆ φX ˆ φY , ˆ φX,Y −ˆ φX ˆ φY  H X H Y =  ˆ φX,Y −φX,Y +φX φY −ˆ φX ˆ φY , ˆ φX,Y −φX,Y + φX φY −ˆ φX ˆ φY  . (10) Then, a natural bootstrap procedure that allows us to estimate the p-value for the independence test is developed as follows: 1. To randomly sample with replacement n elements from the original dataset D , repeating m times. We denote by D j ∗= { (x j ∗ i , y j ∗ i , r j ∗ i ) } n i =1 , j = 1 , . . . , m the j th random sample obtained. 2. To calculate  HSIC j ∗ ( ˆ P X,Y , ˆ P X ˆ P Y ) as  HSIC j ∗ ( ˆ P X,Y , ˆ P X ˆ P Y ) =  ˆ φX,Y −ˆ φj ∗ X,Y + ˆ φj ∗ X ˆ φj ∗ Y −ˆ φX ˆ φY , ˆ φX,Y −ˆ φj ∗ X,Y + ˆ φj ∗ X ˆ φj ∗ Y −ˆ φX ˆ φY  H X H Y , (11) where j = 1 , . . . , m , ˆ φj ∗ X,Y (·) , ˆ φj ∗ X (·) , and ˆ φj ∗ X (·) are the kernel mean embeddings estimated from the j th bootstrap sample D j ∗= { (x j ∗ i , y j ∗ i , r j ∗ i ) } n i =1 . 3. To estimate the p-value as p -value = 1 m m  j=1 I  HSIC j ∗ ( ˆ P X,Y , ˆ P X ˆ P Y ) ≥ HSIC ( ˆ P X,Y , ˆ P X ˆ P Y ) . (12) Bootstrap consistency with missing data can be proved by using standard tools of empirical process theory [16] , and it is provided elsewhere. 2.3. Variable selection Independence screening methods select predictor variables based on individual prediction ability; hence, they are ineffective in selecting a subset of variables that are individually weak but strong in combination. Subset selection aims to overcome this drawback by considering and evaluating the prediction ability of a subset of variables as a whole. One popular approach to subset selection is to directly optimize an objective function consisting of two terms: a data fitting term to attain prediction accuracy and a regularization term to penalize a large number of variables [17] . Subset selection has recently been approached using the RKHS paradigm with satisfactory results. Two strategies stand out: first, minimizing the trace of the conditional covariance operator [18] and second, identifying those variables with a non-zero gradient function [19] . The first strategy scales poorly with the number of variables used. The second strategy can be formulated in a more compact manner. Here, it is extended to missing data. Following [19] , we identify the relevant predictors by learning the gradient of regression function f. Thus, it is assumed that if variable X r is not relevant for predicting Y , then g r = ∂ f(X ) /∂ X r = 0 for any value of X . Let us denote by g (X ) = ∇f(X ) = (g 1 (X ) , . . . , g p (X )) T the gradient function. In a small neighborhood of x i we can use the Taylor expansion to approximate f(X ) , so when x j is sufficiently close to x i , f(x j ) ≈y i + g(x i )(x j −x i ) . We then define the estimation error as a function of g (·) : E(g ) = E X ,Y, X  ,Y  [ ω(X , X  )(Y −Y  −g (X ) T (X −X  ))] 2 , where X  , Y  denote independent and random variables distributed as X and Y , respectively. Function ω(x i , x j ) is an appropriate weight function that decreases as  x i −x j  increases and ensures that the local neighbourhood of x i contributes more to estimating the gradient g (x i ) . Typically, ω(x i , x j ) = e − x i −x j  2 /τ2 n , where τ2 n is a positive parameter which must be adjusted to ensure asymptotic estimation consistency. Because only a limited number of samples { (x i , y i , r i ) } n i =1 are observed, we approximate E(g ) using its empirical counterpart ˆ E (g ) = 1 n 2 n  i,j=1 ω ij (y j −y i −g (x i ) T (x j −x i )) 2 , (13) where ω ij = ω(x i , x j ) . We can add a regularization term for enforcing a sparsity constraint on the gradient vector, with the aim of shrinking the partial derivatives g r towards zero with respect to irrelevant variables. We then add the term J(g ) = λn  p r=1 ηr J(g r ) , where ηr are adaptive tuning parameters. On the other hand, we can define the estimation error in (13) as a functional in the RKHS H p k , and thus g ∈ H p k and E : H k ×p ···×H k → R + , induced by a pre-specified positive kernel k . Thus, we propose the following optimization formula to learn the gradient vector: arg min g ∈H p k 1 n 2 n  i,j=1 ω ij (y i −y j −g (x i ) T (x i −x j )) 2 + J(g ) . (14) 4 M. Matabuena, P. Félix, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 221 (2022) 106905 Under the MAR assumption, we propose substituting ω ij weights with ˆ ω ∗ ij = ˆ ω ∗ i ˆ ω ∗ j ω ij , where ˆ ω ∗ i and ˆ ω ∗ j denote the estimated normalized weights associated with data i th and j th , respectively, according to (8) . The variable selection expression can be rewritten as follows: arg min g ∈H p k 1 n 2 n  i,j=1 ˆ ω ∗ ij (y i −y j −g (x i ) T (x i −x j )) 2 + J(g ) . (15) The representer theorem states that the minimizer of (15) can be represented as a finite linear combination of kernel products evaluated on the dataset samples [20] : g r (·) = n  i =1 αr i k X (·, x i ) , r = 1 , . . . , p, (16) where αr ∈ R n . Given this representation, g r (·) = 0 iff αr = (αr 1 , ... , αr n ) T = (0 , . . . , 0) T , or more concisely,  αr  2 = 0 . Several regularization terms have been considered in previous studies. We adopted the group lasso penalty [19,21] : J(g r ) = inf   αr  2 : g r (·) = n  i =1 αr i k X (·, x i )  , (17) which encourages all αr i , i = 1 , . . . , n to be selected or shrunk to zero together to achieve the purpose of variable selection. Thus, our optimization problem can be rewritten as: arg min α1 , ... ,αp n  i,j=1 ˆ ω ∗ ij (y i −f ∗(x i , x j )) 2 + λn p  r=1 ηr  αr  2 , (18) where f ∗(x i , x j ) = y j − p r=1 k T i αr (x r i −x r j ) , k i = (k (x i , x 1 ) , . . . , k (x i , x n )) T is the i th row of K = (k (x i , x j )) n ×n , and λn are tuning parameters. This last expression simplifies the original optimization framework (14) from a functional space to a vector space, and it can be solved in O (| U| 2 p 2 ) using a block coordinate descent algorithm [19] . 2.4. Prediction and uncertainty analysis Let us recall that the ultimate goal is to predict Y using the information provided by predictor variables X . To achieve this aim, we adopt the kernel ridge regression approach proposed by Liu and Goldberg [22] . However, we draw on the linear regression theory to efficiently compute the leave-one-out cross-validation regularization parameter. This class of regularization parameters has been proven to largely shape the model performance [23] . Furthermore, estimating the uncertainty of predictions by providing robust confidence intervals is a valuable tool for subsequent decisions. Thus, we compute intervals with good finite sample coverage using advances in conformal inference recently exploited in causal theory [24] . Let us assume a linear regression model: y i = f(x i ) + = x T i β+ i = 1 , . . . , n, (19) where βis the vector of coefficients of the linear model. Given the original dataset D = { (x i , y i , r i ) } n i =1 , kernel ridge regression is based on solving the following optimization problem: ˆ β= arg min β∈ R p n  i =1 (y i −x T i β) 2 + λ|| β|| 2 2 , (20) which is solved by ˆ β= (x T x + λI ) −1 x T y , where x = (x 1 , . . . , x n ) T , y = (y 1 , . . . , y n ) T , and λ> 0 is the smoothing parameter of the regularization term. Let H k be an RKHS with a kernel k X . Then, by replacing every x i with φ(x i ) and assuming that β=  n i =1 αi φ(x i ) , we obtain an analogue solution to that of Eq. (20) by exploiting the linear structure of the problem but changing the usual dot product by the inner product of the selected RKHS. In particular, ˆ α= (K + λI ) −1 y , where K = (k (x i , x j )) n ×n . In [22] , the authors proposed two estimators for the missing data. In both cases, the solution has the same closed-form expression, given by the representer theorem. The first is ˆ α= (λI + W ) −1 Wy , where the missing data mechanism is handled using the IPW estimator. The second is obtained through doubly robust estimation, combining a preliminary imputation of the missing response with the IPW estimator: ˆ α= (K + λI ) −1 (Wy + (I −W ) μ(x )) , (21) where W = diag(w 1 , . . . , w n ) denotes a diagonal matrix containing the weights (see Eq. (7) ) and μ(x ) = (μ(x i ) , . . . , μ(x n )) T denotes the imputation function. Doubly robust estimators achieve optimal asymptotic variance when their weights w 1 , . . . , w n and imputation function μ(·) are correctly specified, and only one of them needs to be correctly specified to achieve consistency. However, when any of them fails, the regression model performance can deteriorate dramatically with a finite sample [25,26] , thereby failing to provide real advantages with respect to the IPW estimator. The impact of the smoothing parameter on model generalization is an essential issue for the ensuing performance and is strongly related to the minimum-norm interpolation problem in the context of RKHS. Therefore, we propose the selection of the smoothing parameter through leave-one-out cross-validation by adapting estimators to missing data [23] . To supply a prediction interval for the response with a confidence level of 1 −α, we provide a novel algorithm for performing conformal inference [24,27] , which is valid for handling missing responses and heteroscedastic noise. We randomly split the dataset D = { (x i , y i , r i ) } n i =1 into training and test sets D train = { (x train i , y train i , r train i ) } n 1 i =1 and D test = { (x test i , y test i , r test i ) } n 2 i =1 , where n = n 1 + n 2 . For a given new observation x n +1 we go through the following steps: 1. Fit the mean regression function ˆ f (·) from the set D train , according to Eq. (21) . 2. Compute the residuals ˆ i = | y test i −ˆ f (x test i ) | / ˆ σ(x test i ) , for every i = 1 , . . . , n 2 with r test i = 1 . The value ˆ σ(x test i ) is estimated by a regression function that predicts the absolute deviation of the residuals fitted with the training sample. 3. Estimate the empirical distribution as follows: ˆ F  n 2 +1 (x ) = 1  n 2+1 i =1 ˆ w test i  n 2  i =1 1 { ˆ i ≤x } ˆ w test i + ˆ w test n 2+1  , (22) where we also incorporate the weights of x n +1 and ˆ w test n 2+1 into the estimation. 4. Compute the 1 −αquantile, ˆ q 1 −α, from ˆ F  n 2 +1 . 5. Finally, return [ ˆ f (x n +1 ) −ˆ q 1 −αˆ σ(x n +1 ) , ˆ f (x n +1 ) + ˆ q 1 −αˆ σ(x n +1 )] as the required prediction interval. 3. An application in modelling long-term changes in glucose levels Diabetes mellitus is one of the most critical public health problems and the ninth major cause of mortality worldwide [28] . At present, over 416 million and 47 million patients have type 2 (T2D) and type 1 (T1D) diabetes, respectively, [29] . Significantly, around 50% of patients with diabetes are undiagnosed [29] . Considering the impact of this pandemic among the general population, there is a need for new health policies and guidelines to enable early 5 M. Matabuena, P. Félix, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 221 (2022) 106905 recognition of at-risk patients and improvement in the methodology of disease diagnosis in the standard clinical routine [30] . Some previous studies have focused on developing predictive models for patient stratification. Thus, the Finnish FINDRISC provides a diabetes score to predict the probability of developing diabetes within ten years using logistic regression [31] . In addition, the German GDRS provides a different score to predict the time to becoming a diabetic person using a survival model based on Cox regression [32] . In contrast, some authors argue against using thresholds and categorising patients into different ranges of glucose levels, and hence, against defining diabetes as a homogeneous disease [33] . The availability and rapid adoption of new digital medical devices have enabled an emerging clinical paradigm based on precision medicine, which will be called to raise early diagnosis and guide subsequent clinical decision-making through the intensive use of statistical models and machine learning techniques [34– 37] . In the case of diabetes, the latest advances in sensing technology allow for the assessment of glucose metabolism at a highresolution level by capturing the individual differences in glucose fluctuations at different time scales via continuous glucose monitoring (CGM) [38] . In this sense, although T1D cannot be prevented at present, monitoring is of utmost importance. Recent studies have shown improved glycemic control and decreased rates of hypoglycaemia in T1D patients using CGM, leading both the Endocrine Society and the American Diabetes Association to state that CGM use represents the standard of care for T1D [39,40] . With respect to T2D, strong scientific evidence shows that it can be prevented by regular exercise, healthy eating, and the control of blood pressure and lipids [41] , spurring innovation in wearable technology to enable its prediction and prevention in the general population [42] . Still, few studies have explored the use of CGM data from healthy populations to draw new conclusions regarding glucose homeostasis. It is worth mentioning [43] , which provides some remarkable insights into the heterogeneity of glucose dysregulation, highlighting the inadequacy of a common designation as T2D for categorising and subsequently managing predictably different conditions. Importantly, this study refutes the assumption of similar glucose excursions for the same amount and composition of food. Ultimately, the specific glucose profiles observed for each patient depend on the complex interplay between the pathophysiological mechanisms of insulin resistance and insulin secretion [38,43] , thus enabling the treatment of the glycemic profile of an individual as a personal signature of glucose homeostasis. Accordingly, an appropriate interpretation of CGM data could help identify early stages of glucose dysregulation in apparently healthy individuals, with the possibility of providing early and tailored interventions. In this sense, there is a need for further research on the predictive value of CGM data [38] . This study aimed to examine the predictability of long-term changes in glucose levels by using a random sample of the general population. The exploration of the predictive value of the information provided by CGM data is of particular interest. For this purpose, we consider a new distributional representation and compare the corresponding model performance with well-established biomarkers for the diagnosis and management of diabetes. 3.1. The AEGIS diabetes study The AEGIS diabetes population study, conducted in the Spanish town of A Estrada (Galicia), aimed to analyze the longitudinal changes in some clinical features related to circulating glucose in 1516 patients over 5 years. In addition, non-routine medical tests, such as CGM, are performed every five years on a randomized sub- Table 1 Characteristics of AEGIS study participants with CGM monitoring by sex. Means and standard deviations are shown. A1c: glycated haemoglobin; FPG, fasting plasma glucose; HOMA-IR, homeostasis model assessment-insulin resistance; BMI, body mass index; CONGA, glycemic variability in terms of continuous overall net glycemic action; MAGE, mean amplitude of glycemic excursions; MODD, mean of daily difference. Men (n = 220) Women (n = 361) Age, years 47 . 8 ±14 . 8 48 . 2 ±14 . 5 A1c, % 5 . 6 ±0 . 9 5 . 5 ±0 . 7 FPG, mg/dL 97 ±23 91 ±21 HOMA-IR, mg/dL. μIU/mL 3 . 97 ±5 . 56 2 . 74 ±2 . 47 BMI, kg/m 2 28 . 9 ±4 . 7 27 . 7 ±5 . 3 CONGA, mg/dL 0 . 88 ±0 . 40 0 . 86 ±0 . 36 MAGE, mg/dL 33 . 6 ±22 . 3 31 . 2 ±14 . 6 MODD 0 . 84 ±0 . 58 0 . 77 ±0 . 33 set composed of 581 patients. At the beginning of this study [44] , 581 participants were randomly selected to wear a CGM device for 3 −7 days. Of the 581 participants, 68 were diagnosed with diabetes before the start of the study and 22 were diagnosed during the study. Table 1 lists the baseline characteristics of the 581 patients grouped by sex. After a five-year follow-up, a significant fraction of those individuals did not agree to perform a second glucose monitoring, while some five-year relevant outcomes such as A1c could only be measured in 339 patients. Complete details of the study design and measurement methodology protocol can be found in [44] . 3.2. A distributional representation for CGM data: Glucodensity We adopted a novel functional representation for CGM data, termed glucodensity, which allows us to obtain a personalized functional profile of patient glucose homeostasis [10] . Glucodensity is a natural extension of those metrics based on the time spent in certain ranges. Time in range (TIR) measures the proportion of time a person spends with their blood glucose levels within the target range of 70 −180 mg/dL. Accordingly, the time below range (TBR) and time above range (TAR) are defined. These are the current gold standard for representing CGM data [45,46] . Although very intuitive, they have two main disadvantages: first, the range fits poorly depending on the characteristics of the population examined; second, there is a loss of information caused by the discretization of the recorded data into intervals. Instead, glucodensity effectively measures the proportion of time each individual spends at a specific glucose concentration. Given a series of CGM data { x j } m j=1 , the glucodensity can be modeled as a probability density function f(·) that can be approached by kernel density estimation: ˆ f (x ) = 1 m m  j=1 1 h k x −x j h , (23) where h > 0 is the smoothing parameter and k (·) denotes a nonnegative real-valued integrable function ( Fig. 2 ). A relevant issue in kernel analysis is to measure the difference between two density functions. In this study, we use the 2- Wasserstein distance. Given two glucodensities, ˆ f 1 and ˆ f 2 , the 2- Wasserstein distance between them is given by d W 2 ( ˆ f 1 , ˆ f 2 ) =   1 0 ˆ Q f 1 (t) −ˆ Q f 2 (t)  2 dt, (24) where ˆ Q f 1 , and ˆ Q f 2 are the corresponding quantile functions. Because the 2-Wassertein distance between two densities depends only on their quantile functions, it is not necessary to resort to 6 M. Matabuena, P. Félix, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 221 (2022) 106905 Fig. 2. Glucodensities estimated from a random sample of the AEGIS study on diabetic and normoglycemic patients. Left: Representation of the proportion of time spent by a patient at each glucose concentration over a continuum is shown. Right: Representation of the glucodensities in the space of quantile functions is shown. A meaningful difference between the average profiles for both groups of patients can be noticed. density estimation methods, and we can approximate this distance using quantile-function estimations through empirical distributions. Intuitively, glucodensity is more sensitive than the previous CGM summary metrics. We then explored its use in modeling long-term glucose changes and compared it with TBR, TIR, and TAR measures. In addition, we also considered different summary metrics derived from CGM data [47,48] : CONGA (continuous overall net glycemic action), MAGE (mean amplitude of glycemic excursions), and MODD (mean of daily differences). 3.3. Integrating multiple data sources RKHS offers a powerful and natural data analysis paradigm that can cope with data of different natures [49] . A crucial issue is to select a suitable kernel that accurately captures the differences and specific characteristics of each information source examined. In our particular case, we take into account a continuous probability distribution and certain real-valued and categorical data x = (x gluco , x real , x categ ) . A reasonable choice commonly used in the literature is the Laplacian kernel, K = (k (x i , x j )) n ×n , where k (x i , x j ) = e −|| x i −x j || σ. Here, we propose using the Laplacian kernel with the standard Euclidean distance as a universal and characteristic kernel in a real vector space. Moreover, it can be shown that the Laplacian kernel retains these properties considering the set of continuous density functions endowed with 2 −Wasserstein distance, providing theoretical guarantees that we can approximate a large variety of regression functions. Based on the connection between positive kernels and negative-type metrics [50,51] , we propose using a simple and global Laplacian kernel that integrates these three sources: k X (x i , x j ) = e −a || x gluco i −x gluco j || σgluco + b || x real i −x real j || σreal + c || x categ i −x categ j || σcateg , (25) where a, b, c, σgluco , σreal , σcateg > 0 and we assume for the sake of simplicity that (a, b, c) ∈ S 2 , where S 2 is a 2-simplex, S 2 = { (a, b, c) ∈ R 3 : a + b + c = 1 ;0 ≤a ≤1 , 0 ≤b ≤1 , 0 ≤c ≤1 } . 4. Results The present framework of predictive tools allows us to answer some open clinical questions concerning long-term glucose changes from the analysis of data in the AEGIS study. 1. Glycated haemoglobin A1c is a haemoglobin-glucose combination formed within the cell; it is a useful indicator of long-term blood glucose control and is considered the standard biomarker for diabetes diagnosis and management. Is there a prognostic variable that can be used to predict future A1c changes in healthy individuals? 2. Current medical literature assigns a considerable relevance to all of the predictor variables listed in Table 1 for characterizing the evolution and impact of glucose homeostasis on health. However, from a biological perspective, these variables are well known to be highly correlated. Can we identify a reduced subset of relevant explanatory variables to predict five-year A1c changes? 7 M. Matabuena, P. Félix, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 221 (2022) 106905 Table 2 Estimated raw p-values of A1c total variation vs each biomarker using the method proposed in Section 2.2 with normoglycemic patients. Variable p −v alue Age 0.32 Sex 0.16 FPG 0.50 HOMA-IR 0.52 BMI 0.42 A1c 0.03 CONGA 0.24 MAGE 0.68 MODD 0.16 Glucodensity < 0 . 001 3. CGM technology may provide a more suitable tool for assessing glucose homeostasis than traditional diabetes biomarkers. How does CGM data improve our ability to predict future A1c changes? 4. An increased uncertainty in predictions for a specific region of the feature space may suggest a subpopulation that has not been properly modeled. Can we provide a characterization of individuals for whom the model yields a less accurate prediction? 4.1. Is there a prognostic variable that can be used to predict future A1c changes in healthy individuals? To answer this question, we studied whether there was any evidence of a univariate statistical association for normoglycemic patients (A1c < 5.7% and FPG < 100 mg/dL) between glucose variation measured by A1c 5 years −A1c initial and the predictor variables shown in Table 1 . For this purpose, we use the Hilbert-Schmidt independence criterion proposed in the context of missing data ( Section 2.2 ) together with a specific bootstrap approach designed for this task. The underlying mechanism of missing data was estimated using univariate logistic regression. The results in Table 2 show that the only statistically significant variables with a p-value of less than 5% are glucodensity and basal A1c. Fig. 3 illustrates that the marginal relationships with other variables, if any, are weak. 4.2. Can we identify a reduced subset of relevant variables to predict five-year A1c changes? Multivariate models can exploit higher-order interactions between the predictors and the response to improve predictions. However, a key point in increasing the interpretability and generalization ability of the model is to identify a subset of the variables that capture the essential information in the dataset, thus removing redundancy. We adjusted the method proposed in Section 2.3 to find the subset of variables most strongly associated with A1c 5 years . For this purpose, both diabetic and non-diabetic patients were analysed, and we considered all the variables in Table 1 except for sex. We also included the TBR, TIR, and TAR measures specified in Section 3.2 . To avoid overfitting and improve the reproducibility of the results, we selected model parameters by cross-validation. We estimated the underlying missing data mechanism using lasso logistic regression. Finally, the explanatory variables selected by the algorithm were age, A1c initial , FPG, BMI, and MAGE. Notably, the CGM contribution is made through the specific MAGE index, leaving aside time in ranges. 4.3. How does CGM data improve our ability to predict future A1c changes? To answer this question, we fit several kernel ridge regression models ( Section 2.4 ) for predicting A1c 5 years : 1) excluding CGM data as a predictor; 2) including CGM data through the MAGE index; 3) including CGM data through the above-mentioned time in ranges; and 4) including CGM data through glucodensity representation. Both share age, A1c initial , FPG, and BMI as covariates. The kernel selection and parameter tuning were calibrated as described in Section 3.3 . To compare the performance of these regression models, we used R 2 after including the specific missing data mechanism: R 2 = 1 − n i =1 y i − n j=1 ,j = i w j y j  n j=1 ,j = i w j 2  n i =1 y i −ˆ f −i (x i ) 2 , (26) where ˆ f −i (·) , is the regression function fitted to { (x j , y j , r j ) } n j = i , i.e. excluding the i th -datum. The performance results, obtained using leave-one-out crossvalidation, are: 1) R 2 noCGM = 0.61; 2) R 2 MAGE = 0.65, 3) R 2 T IR = 0.64; and 4) R 2 gluco = 0.71. For comparison purposes, partial least squares regression has also been applied, by using leave-one-out crossvalidation, obtaining the following results: 1) R 2 noCGM = 0.52; 2) R 2 MAGE = 0.57; and 3) R 2 T IR = 0.56. Fig. 4 depicts the residuals versus the A1c initial values for the model including glucodensities. As can be seen, the highest residuals are found in diabetic patients; besides that, their distribution is heterogeneous. As expected, the prediabetic range is the most challenging to model. From a total of 60 prediabetic (according to A1c criteria: 5.7% ≥A1c ≥6.4%) individuals at baseline, 22 of them were diagnosed with diabetes (A1c ≥6.5% or FPG ≥126 mg/dL) 2 five years later, 14 of them with A1c ≥6.5%; 5 of the latter were correctly predicted by our models and 9 were wrongly predicted. Because of the limited number of individuals in this range any extrapolation needs to be restrained. Ultimately, CGM data represented by glucodensities provide valuable information for predicting long-term A1c changes. 4.4. Can we provide a characterization of individuals for whom the model yields a less accurate prediction? Fig. 5 depicts prediction intervals at a confidence level of 90 % , after applying conformal inference ( Section 2.4 ) to measure the uncertainty of the predictions performed by the above regression model (CGM data included as a covariate). We regard an A1c 5 year prediction as significantly affected by uncertainty if the length of the interval is greater than 0.7 because a deviation greater than this threshold can entail a change in the glycemic state of the patient, for example, from normoglycemic to diabetes. Hence, we can identify certain clinical features that allow us to assign each patient to high-or low-variability groups based on the uncertainty of future glucose values. This can be useful to phenotypically characterize some subpopulations for which the model provides an unreliable prediction, and therefore, a more personalised follow-up is advisable. In particular, Fig. 6 shows that long-term changes cannot be adequately predicted for individuals with elevated FPG levels. The same holds true for individuals with FPG levels in the normoglycemic range and overweight. More refined decision rules can be established at higher measurement costs. 2 Diagnostic criteria according to American Diabetes Association. 8 M. Matabuena, P. Félix, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 221 (2022) 106905 Fig. 3. Marginal dependence relation between examined variables in the AEGIS database. 5. Discussion The above analysis aimed to explore the predictability of glucose regulation in the general population by studying the relationship between patient basal characteristics at the start of a longitudinal study and A1c values obtained five years later. Specifically, we intend to exploit the ability of CGM data to effectively capture a personal signature of glucose homeostasis through the inclusion of glucodensity, a novel distributional representation of glucose excursions, as a predictor. The AEGIS study makes it possible to assess the predictive capacity of glucodensity in the context of well-known biomarkers for diabetes diagnosis and control, providing some interesting findings. First, glucodensity shows a significant association with A1c changes, using statistical dependence measures in normoglycemic subjects. Nevertheless, the weak marginal association of biomarkers with A1c 5 years suggests the need for a multivariate approach to capture the complexity of long-term glucose changes. The application of a variable selection procedure supplies a subset of relevant biomarkers (age, A1c initial , FPG, BMI, and MAGE) resulting from the detection of higher-order interactions with A1c 5 years . Then, the ability to predict A1c 5 years from this subset of biomarkers is analysed using several regression models that differ in terms of including CGM data as a predictor. As a result, the R 2 gluco value, corresponding to the model which adopts a glucodensity-based representation for CGM data, shows a good proportion of variance explained by the model and is similar to that reported by other authors for short-term predictions [52,53] . Moreover, glucodensity has a positive impact on improving accuracy in predicting A1c 5 years by expanding the model expressiveness along the continuous spectrum of glucose concentrations. Some recent studies have proposed different machine learning methods for predicting the progression to diabetes from a healthy or prediabetic state with relatively good performance [54,55] . The strength of both studies lies in their inclusion of a large number of subjects. Both of them proposed a classification strategy relying on a threshold-based categorization upon different ranges of glucose or A1c. This allows them to obtain robust results but at the expense of making an inappropriate interpretation in physiological terms [33,38] . Furthermore, both studies are of observational nature and subjects with only partially available data were excluded from the analysis. This suggests that caution should be exercised when utilising these results for decision-making. The present line of research assigns a key role to the analysis of glucose excursions from CGM data in search of better phenotyping and corresponding progress towards the implementation 9