Modeling conditional reference regions: application to glycemic markers

Received: 22 June 2020 Revised: 7 July 2021 Accepted: 29 July 2021 DOI: 10.1002/sim.9163 RESEARCH ARTICLE Modeling conditional reference regions: Application to glycemic markers Óscar Lado-Baleato1Javier Roca-Pardiñas2Carmen Cadarso-Suárez1 Francisco Gude3 1Department of Statistics, Mathematical Analysis, and Optimization, Universidade de Santiago de Compostela, Galicia, Spain 2Statistical Inference, Decision and Operations Research, Universidade de Vigo, Galicia, Spain 3Clinical Epidemiology Unit, Complexo Hospitalario de Santiago de Compostela, Galicia, Spain Correspondence Óscar Lado-Baleato, Department of Statistics, Mathematical Analysis, and Optimization, Universidade de Santiago de Compostela, Galicia, Spain. Email: [email protected] Funding information Consellería de Cultura, Educación e Ordenación Universitaria, Xunta de Galicia, Grant/Award Numbers: ED341D-R2016/032, ED431C2016-025, IN607B 2018-1; Instituto de Salud Carlos III, Grant/Award Numbers: PI16/01404, RD16/0017/0018; Ministerio de Economía y Competitividad, Grant/Award Number: MTM2017-83513-R Many clinical decisions are taken based on the results of continuous diagnostic tests. Usually, only the results of one single test is taken into consideration, the interpretation of which requires a reference range for the healthy population. However, the use of two different tests, can be necessary in the diagnosis of certain diseases. This obliges a bivariate reference region be available for their interpretation. It should also be remembered that reference regions may depend on patient variables (eg, age and sex) independent of the suspected disease. However, few proposals have been made regarding the statistical modeling of such reference regions, and those put forward have always assumed a Gaussian distribution, which can be rather restrictive. The present work describes a new statistical method that allows such reference regions to be estimated with no insistence on the results being normally distributed. The proposed method is based on a bivariate location-scale model that provides probabilistic regions covering a specific percentage of the bivariate data, dependent on certain covariates. The reference region is estimated nonparametrically and the nonlinear effects of continuous covariates via polynomial kernel smoothers in additive models. The bivariate model is estimated using a backfitting algorithm, and the optimal smoothing parameters of the kernel smoothers selected by cross-validation. The model performed satisfactorily in simulation studies under the assumption of non-Gaussian conditions. Finally, the proposed methodology was found to be useful in estimating a reference region for two continuous diagnostic tests for diabetes (fasting plasma glucose and glycated hemoglobin), taking into account the age of the patient. KEYWORDS conditional reference regions, diabetes, flexible additive predictors, kernel smoothing, regression 1INTRODUCTION The diagnosis and treatment of disease commonly rests on the results of measureable biomarker-based clinical laboratory tests. Indeed, some 70% of the decisions made in clinical practice are taken based on such results.1For every test result, the clinical laboratory provides comparator values to help the clinician understand in context the information provided. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distributioninanymedium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. © 2021 The Authors. Statistics in Medicine published by John Wiley & Sons Ltd. 5926 wileyonlinelibrary.com/journal/sim Statistics in Medicine. 2021;40:5926–5946. LADO-BALEATO et al. 5927 These comparator values are often referred to as the reference interval,2they usually reflect the range of values within which 95% of the results of the normal healthy population falls. When a single biomarker is examined, the reference interval is classically obtained using quantile estimation techniques,3or by conditional quantile regression if any variable modifies the distribution of the response variable (ie, the reference curve).4,5 However, there is often more than one test for diagnosing a disease. For instance, the diagnosis of diabetes may be based on plasma glucose criteria, such as fasting plasma glucose (FPG) or the 2-hour plasma glucose value obtained during a 75-g oral glucose tolerance test, or the glycated hemoglobin (HbA1c) test.6The same tests may be used to screen for, diagnose and monitor the effects of treatment for diabetes. However, measuring glycemic control is not foolproof; its clinical usefulness is affected by a number of biological and analytical factors. Disagreement between glycemic control measurements are common, and clinicians need to know what might explain them.7 Following the guidelines of the American Diabetes Association, the diagnosis of diabetes is defined as FPG levels of ≥126 mg/dL, and of HbA1c ≥6.5%. The same cut-offs are used for children, adolescents, and adults. FPG and HbA1c levels are reported to be strongly correlated,8both in members of the general population and in patients with diabetes.9,10 If two tests are indifferently used for the diagnosis of a disease, the correlation between them should be strong. Therefore, when diagnosing a disease using two markers, it might be reasonable to estimate their combined multivariate reference region instead of the reference interval for each test. The idea of combining reference regions for two or more laboratory tests has been discussed in the biomedical,11-13 and statistical literature.14 However, the proposals made have required responses that follow a multivariate Gaussian distribution, condition that are not fulfilled by many markers used in the clinical setting. For instance, FPG and HbA1c concentrations both show a skewed distribution. Hence, a more general method for estimating multivariate reference regions is needed. Moreover, since HbA1c increases with age even after adjusting for glucose levels,15-17 this variable should be taken into account when establishing cut-offs or reference values. The literature reports but a few attempts to define reference regions for non-Gaussian multivariate responses. Nonparametric reference regions for multivariate responses can be estimated using multivariate quantiles. However, there is no single definition of what a multivariate quantile is.18 Most current definitions are based on a center-outward ordering of the data points,19,20 defining a convex hull in which a proportion of the more central data points falls. Halfspace depth bivariate quantiles based on directional projections have recently been extended to the regression setting for estimating bivariate contours conditioned by covariates.21 These conditional bivariate quantile models have been used in the study of anthropometric characteristics affected by age.22 However, the clinical interpretation of the results is not clear. Easily interpretable conditional bivariate quantile contours can be estimated using the Wei model,23 which is based on the estimation of i) the marginal stratified conditional quantile regression for each response, and ii) bivariate quantiles using simulated data points from the above marginal stratified conditional quantile regression. The main drawback of this proposal is the influence of the univariate quantile regression outcome on the performance of the final bivariate quantile contour. Another nonparametric conditional method for detecting extreme combinations of two variables has also been proposed.24 This does not define a reference region, but four bivariate reference curves for detecting four possible atypical combinations of two variables. However, this alternative returns a higher false discovery rate for the reference region. Non-Gaussian bivariate regions for detecting joint outliers may be estimated using conditional copula regression models.25 This procedure was recently applied to identify children with abnormal vision.26 Copula regression models allow a bivariate distribution to be constructed from a modular perspective, combining two univariate parametric distributions and the parametric copula joining them.27 Finally, the response parameters are made dependent on the covariates using flexible additive models.28 Given copula regression model structure, and the effect of the estimated covariate on the means, variances and correlation of the responses, a conditional bivariate region for exploring the atypical combinations of two measurements can be estimated. The parametric representation of these models means several choices have to be made in the model building process, which complicates the use of this methodology when dealing with real data problems. The present work proposes a new regression model for estimating conditional bivariate reference regions. The reference region is defined as the density function contour level which contains the bivariate data points with a given probability depending on the covariates. Unlike existing methods, cited above, our statistical model places no parametric restrictions on the response, and the nonlinear effects of continuous covariates may be estimated using local polynomial regression smoothers. Our proposal is an extension to bivariate data of a previous work,29 where the authors used a location-scale model to estimate univariate percentile curves. The final performance of our conditional reference regions depends heavily on a bivariate kernel density estimator. It is well known, that nonparametric density estimators are largely influenced by kernel bandwidth matrix choice. Therefore, in this work, we also propose a new kernel bandwidth estimation method in order to obtain a reference region with a coverage of the bivariate data points close to the nominal level. 5928 LADO-BALEATO et al. The proposed statistical methodology was tested with simulated data, and applied to estimate an age dependent reference region for two glycemic markers (FPG and HbA1c), in a cohort of nominally normoglycemic persons. The remainder of this article is organized as follows. Section 2 introduces the model structure and its estimation details. Section 3 evaluates the conditional reference region estimation, contemplating scenarios in which the model could be partially miss-specified. Section 4 makes use of the model for estimating an age-specific bivariate reference region for the diabetes markers FPG and HbA1c. Finally, Section 5 discusses the limitations and applicability of the proposed model, highlighting possible extensions. 2CONDITIONAL BIVARIATE REFERENCE REGION MODEL This section presents a regression model for estimating a bivariate reference region conditioned by a set of covariates. The proposed model has no parametric restriction, and is based on the estimate of a location-scale regression, taking into account the effect of the covariate on the correlation between the response variables. 2.1 Model formulation Let X=(X1,…,Xp)be a vector of p covariates, and let Y=(Y1,Y2)be a continuous bivariate response of interest. In this context, the aim is to obtain a bivariate region 𝜏of Yconditioned by the covariates Xdenoted as R𝜏(X), and containing the 100𝜏%of the bivariate data points. The following structure is assumed: Y=(𝜇1(X) 𝜇2(X))+𝚺1∕2(X)(𝜀1 𝜀2),(1) where 𝜇1and 𝜇2represent the response means, and the variance-covariance matrix is given by: 𝚺(X)=(𝜎2 1(X)𝜎12(X) 𝜎12(X)𝜎2 2(X)) The bivariate residuals (𝜀1,𝜀 2)are assumed to be independent of the covariates and to have a mean of zero, zero unit variance, zero correlation, and an unknown density function f(𝜀1,𝜀 2). Note that, Σ1∕2(X)represents the Cholesky decomposition of the variance-covariance matrix Σ(X)so that, Σ1∕2(X)(Σ1∕2(X))T=Σ(X). Thus, for any given Xthe bivariate region for (Y1,Y2)is given by: R𝜏(X)=(𝜇1(X) 𝜇2(X))+𝚺1∕2(X)R𝜏for 𝜏∈[0,1],(2) where the R𝜏is the unconditional bivariate region containing the 100𝜏%of the model residuals (𝜀1,𝜀 2)defined as: R𝜏={(u,v)∈R2|f(u,v)≤k}for 𝜏∈[0,1], where kvalue is chosen so that P((𝜀1,𝜀 2)∈R𝜏)=𝜏for 𝜏∈[0,1]. In Equation (1), the response parameters (𝜇1,𝜇 2,𝜎2 1,𝜎2 2,𝜎 12)are related to the covariates vector Xvia additive predictors and known link functions G, which ensure that the restrictions on the parameter spaces are maintained. The following additive predictors are considered: 𝜇r(X)=𝛼r+ p ∑ j=1 fjr(Xj)and 𝜎2 r(X)=G𝜎(𝛽r+ p ∑ j=1 gjr(Xj))for r=1,2,(3) where fjr and gjr are smooth unknown functions, 𝛼rand 𝛽rare intercept coefficients and the link function G𝜎=exp(⋅)to ensure that 𝜎2 r(X)≥0. Finally, the following additive structure is assumed for the responses’ association with one another, LADO-BALEATO et al. 5929 expressed as a linear correlation coefficient (𝜌): 𝜌(X)=G𝜌(𝛾+ p ∑ j=1 mj(Xj)),(4) where mjare smooth unknown functions, 𝛾is an intercept coefficient, and in this case the link function G𝜌=tanh(⋅)to ensure that 𝜌(X)∈[−1,1]. For the sake of mathematical notational simplicity, only a nonlinear effect of the continuous covariates is contemplated in Equations (3) and (4), but they could easily be adapted to incorporate factor effects. For instance, if the first p1covariates define categories, the expression of 𝜇r(X)in (3) can be replaced by the following semiparametric structure 𝜇r(X)=𝛼r+∑p1 j=1𝛼jrXj+∑p j=p1+1fjr(Xj). The variances and correlation structures can be similarly treated. 2.2 Estimation algorithm This section discusses the procedure for estimating the conditional bivariate region presented in Equation (2). The methodology is based on the estimate of the covariate effects on the means of the responses via the use of a flexible additive predictor, and then on the variance-covariance matrix using the squared residuals of the conditional means estimates. Finally, the bivariate region 𝜏is obtained using a bivariate kernel estimate of the standardized bivariate residuals density. Specifically, given a sample of {(Yi1,Yi2),Xi}n i=1,whereXi=(Xi1,…,Xip), the proposed estimation algorithm is as follows: Step 1: For r=1,2 additive predictor is fitted to the original sample {Yir,Xi}n i=1to obtain the estimates: 𝜇r(Xi)=𝛼r+ p ∑ j=1  fjr(Xij)for i=1,…,n.(5) Step 2: For r=1,2 the squared residuals of the previous models 𝜇r(Xi)are obtained and additive predictor fitted to the sample {(Yir −𝜇r(Xi))2,Xi}n i=1and obtain the estimates: 𝜎2 r(Xi)=G𝜎( 𝛽r+ p ∑ j=1  gjr (Xij))for i=1,…,n.(6) Step 3: For i=1,…,nthe standardized residuals are computed:  ri=(Yi1−𝜇1(Xi))( Yi2−𝜇2(Xi)) 𝜎1(Xi)𝜎2(Xi) and the correlation model 𝜌(X)using the sample { ri,Xi}n i=1is fitted as: 𝜌(Xi)=G𝜌(𝛾 + p ∑ j=1  mj(Xij)).(7) Step 4: For i=1,…,nthe estimated standardized bivariate residuals are then obtained: ⎛⎜⎜⎜⎝ 𝜀i1 𝜀i2⎞⎟⎟⎟⎠ = Σ−1∕2(X)(Yi1−𝜇1(Xi) Yi2−𝜇2(Xi)), where  Σ−1∕2(X)is the inverse of the Cholesky decomposition of Σ(X). Using the sample {( 𝜀i1,𝜀i2)}n i=1the kernel estimation of the bivariate density  f(𝜀1,𝜀 2)is obtained. Then, the bivariate region on the residual scale (R𝜏)canbe obtained as 5930 LADO-BALEATO et al.  R𝜏={(u,v)∈R2| f(u,v)≤ k},(8) where  kis the quantile 𝜏of  f(𝜀1,𝜀2). Step 5: Finally, the conditional bivariate region for each Xivalue is given by:  R𝜏(Xi)=(𝜇1(Xi) 𝜇2(Xi))+ 𝚺−1∕2(Xi) R𝜏. Full details about nonlinear functions ( fjr, gjr, mj) estimation can be found in Appendixes A and B. Moreover, in Appendix D, we explain how to obtain confidence intervals for the additive predictors covariates effects through bootstrap. In addition to the R𝜏(X)estimate, the proposed algorithm can be used to obtain the univariate conditional reference curves for each response variable, applying the following expression:  Q𝜏r(X)= 𝜇r(X)+ 𝜎r(X)𝜀𝜏rfor r=1,2,(9) where 𝜀𝜏ris the empirical 𝜏-quantile of the univariate errors 𝜀1r,…,𝜀 nr. 2.3 Bivariate kernel estimation: Bandwidth selection problem The covariate dependent bivariate reference region R𝜏(Xi)requires the estimation of a region containing the standardized bivariate residuals with a given probability 𝜏derived from the bivariate residuals’ density (see Equations (11) and (12)). Given the sample {(𝜀i1,𝜀i2)}n i=1,thefdensity estimator at a given point (u,v)is given by:  f((u,v),H)=1 n n ∑ i=1 KH(u−𝜀i1 v−𝜀i2), where K(⋅)represents the kernel (a bivariate symmetric probability density function, usually the standard bivariate Gaussian distribution), and Ha diagonal matrix defining the kernel bandwidth, the selection of which is crucial for obtaining good estimate of R𝜏and hence the final region R𝜏(Xi). For the selection of the bandwidth H, a plug-in30 or cross-validation31 estimator can be used, as in any density estimation problem. However, the optimal bandwidth for density estimation might not be optimal for the coverage properties of the conditional bivariate region  R𝜏(Xi)(see Figure 1 and Appendix C). Thus, a bandwidth needs to be chosen such that it minimizes the difference between the estimated and nominal coverage of  R𝜏(Xi).Thisisachievedherebybasingthe selection on the expression  H= h𝜆,where his the former bandwidth estimate obtained using the plug-in estimator, and 𝜆is a parameter modulating the final shape of the estimated region. This parameter is estimated as  𝜆=arg min 𝜆||||||(n−1 n ∑ i=1 I{(Yi1,Yi2)∈R(−i)(Xi)})−𝜏|||||| ,(10) where 𝜏is the desired coverage and  R(−i) 𝜏(Xi)is the estimated bivariate region without the ith observation. Given the high-computational cost of (10), a k-fold cross-validation scheme could be used instead. 3SIMULATION STUDY In this section, the estimated conditional bivariate regions  R𝜏(X)are evaluated in terms of the data point coverage and its proximity to the theoretical R𝜏(X). This coverage is estimated using an out-sample design, using one dataset in the estimation and the other in the evaluation. The comparison between the estimated and theoretical bivariate regions was performed using the root mean square error (RMSE) distance. Taking into account that both the theoretical and the estimated bivariate reference regions are defined by a set of bivariate points, R𝜏(Xi)={(ui k,vi k)}B k=1,and R𝜏(Xi)={(ûi j, vi j)}b j=1, LADO-BALEATO et al. 5931 −20246 −20246 Y1|x=0.5 Y2|x=0.5 Theoretical Estimated −20246 −20246 Y1|x=0.5 Y2|x=0.5 Theoretical Estimated FIGURE 1 Example estimating the R𝜏(Xi)for 𝜏=0.10,0.50 and 0.90, using the plug-in bandwidth estimator (left) and following the method proposed in Equation (10) (right). The theoretical bivariate regions were obtained using the parametric density function of the response variable, and 100 000 simulated data points [Colour figure can be viewed at wileyonlinelibrary.com] the RMSE distance was defined as follows: RMSE ( R𝜏(Xi),R𝜏(Xi))=b−1 b ∑ j=1 min 1≤k≤B√(ûi j−ui k)2+(  vi j−vi k)2.(11) In practical applications, the bivariate kernel contour lines are used to estimate (ûi j, vi j). 3.1 Scenario 1: Response with Gaussian error In the first simulation scenario, the datasets were generated according to Equation (1). The means, the variances, and correlation structures of the responses were made dependent on a continuous and a binary regressor. Given a continuous X1∈U[0,1], and a binary covariate X2, the following predictor structure was considered: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝜇1(X)=1+X2+f1(X1) 𝜎2 1(X)=1+X2 𝜇2(X)=1+0.5X2+f2(X1) 𝜎2 2(X)=1+0.5X2 𝜌(X)=0.3+0.2X1+0.3X2 ,(12) where f1(X1)=X1sin(3X1)and f2(X1)=sin(2𝜋X1)represent nonlinear effects. Finally the error term (𝜀1,𝜀 2)was simulated from a standard bivariate Gaussian distribution. The sample size were set to n=500,1000,2000, and the evaluation was done using 1000 replicates. In Figure 2, the R0.95 bivariate region estimates are depicted along with the theoretical estimates for X1=0.5and X2=0. The estimated bivariate region shows a similar shape to the theoretical situation, becoming even closer to it and showing smaller variability as the sample size increases. In Figure 3, the RMSE for  R0.95 is presented for a sequence of continuous predictor values with the binary covariate fixed at zero, for varying sample sizes. As expected, the median RMSE value and variability decreased with increasing sample size. Moreover, the RMSE was higher for the X1values located at the covariate range boundaries (0 and 1); this is justified by the frontier effect associated with nonlinear regression models. 5932 LADO-BALEATO et al. −2 0 2 4 −2024 Scenario 1: n=500 (τ=0.95) Y1 Y2 −2 0 2 4 −2024 Scenario 1: n=1000 (τ=0.95) Y1 Y2 −2 0 2 4 −2024 Scenario 1: n=2000 (τ=0.95) Y1 Y2 −2 0 2 4 −2 0 2 4 Scenario 2: n=500 (τ=0.95) Y1 Y2 −2 0 2 4 −2 0 2 4 Scenario 2: n=1000 (τ=0.95) Y1 Y2 −2 0 2 4 −2 0 2 4 Scenario 2: n=2000 (τ=0.95) Y1 Y2 −202468 −4 −2 0 2 4 6 Scenario 3: n=500 (τ=0.95) Y1 Y2 −2 0 2 4 6 8 −4 −2 0 2 4 6 Scenario 3: n=1000 (τ=0.95) Y1 Y2 −2 0 2 4 6 8 −4 −2 0 2 4 6 Scenario 3: n=2000 (τ=0.95) Y1 Y2 FIGURE 2 Estimated bivariate regions for 𝜏=0.95, X1=0.5andX2=0 along with the theoretical ones (red), for every sample size considered (500, 1000, and 2000) and three simulation scenarios [Colour figure can be viewed at wileyonlinelibrary.com] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 x1=0.0 x1=0.1 x1=0.2 x1=0.3 x1=0.4 x1=0.5 x1=0.6 x1=0.7 x1=0.8 x1=0.9 x1=1.0 RMSE Sample size 500 1000 2000 R0.95(X1) RMSE. Scenario 1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 x1=0.0 x1=0.1 x1=0.2 x1=0.3 x1=0.4 x1=0.5 x1=0.6 x1=0.7 x1=0.8 x1=0.9 x1=1.0 RMSE Sample size 500 1000 2000 R0.95(X1) RMSE. Scenario 2 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 x1=0.0 x1=0.1 x1=0.2 x1=0.3 x1=0.4 x1=0.5 x1=0.6 x1=0.7 x1=0.8 x1=0.9 x1=1.0 RMSE Sample size 500 1000 2000 R0.95(X1) RMSE. Scenario 3 FIGURE 3 Root mean square error of the estimated bivariate regions for 𝜏=0.95, for different sample sizes (500, 1000, and 2000), and with X2=0andX1=(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)for the three simulation scenarios [Colour figure can be viewed at wileyonlinelibrary.com] LADO-BALEATO et al. 5933 Table 1 shows the percentage of bivariate data points contained within the estimated bivariate region to be similar to that seen for the theoretical situation, whether for different nominal levels (𝜏=0.05,0.50,0.90,0.95), continuous predictor variable values, or sample sizes. As for the RMSE, the coverage probability approaches the nominal level as the sample size increases, and becomes worse for values of X1located at the covariate range boundaries (0 and 1). 3.2 Scenarios 2 and 3: Non-Gaussian dependencies In the previous simulation, the proposed model obtained good estimates of the conditional bivariate region. However, one possible problem lies in that, besides the lack of parametric assumptions, the correlation estimate is based on the linear correlation coefficient which is not optimal under some circumstances (non-Gaussian margins and/or nonelliptical structures of dependence). The model was therefore used in two more scenarios contemplating nonstandard bivariate distributions generated from a parametric copula representation. This representation allows complex bivariate distributions to be simulated since the different parametric copula functions proposed in the literature allow for different types of dependence structure between the responses, and each variable may follow any parametric distribution. Given two continuous variables (Y1,Y2), their joint distribution function F1,2conditioned on Xmay be represented in terms of their univariate distributions functions (F1,F2) and a copula Cjoining both, as F1,2(Y1,Y2|X)=C(F1(Y2|𝜇1(X),𝜎 2(X)),F2(Y2|𝜇1(X),𝜎 2(X));𝜃(X)),(13) where F1,F2are two univariate parametric marginal distributions, defined by a location (𝜇) and the scale (𝜎) parameter and Crepresents a parametric copula function defined by an association parameter 𝜃measuring the correlation between both responses. In the present case, the following bivariate response structures were contemplated: •Scenario 2: F1and F2were considered to be two Gaussian distributions and Ca Gumbel copula, defining a nonelliptical structure of dependence (upper-tail dependence). •Scenario 3: F1was simulated from reverse Gumbel distribution and F2from a logistic one.32 As in the former scenario, both are joined by a Gumbel copula. The predictor structures for the bivariate distribution parameters were the same as presented in Equation (12), but the association was expressed in terms of Kendall’s correlation coefficient. The datasets of both scenarios were generated using the gamlss33 and copula34 packages in R. Figure 3 shows the RMSE of the estimated bivariate reference region for both scenarios to decrease with increasing sample size. However, the estimate error is clearly higher than in the first scenario. Nonetheless, the proposed model provides a good approximation to the shape of the theoretical region (see Figure 2), and the percentage of bivariate data points contained within the estimated region is close to the nominal level in both scenarios, becoming better as the sample size increases (see Table 1). These results suggest that the proposed model is quite robust in the face of possible miss-specifications. 4AGE-SPECIFIC BIVARIATE REFERENCE REGION FOR THE GLYCEMIC TESTS 4.1 Motivating database The A-Estrada Glycation and Inflammation Study (AEGIS) is a cross-sectional, population-based study that was performed in the municipality of A Estrada (Galicia, NW Spain). The study objective was to investigate the association between glycation, inflammation status, lifestyle, and common diseases, and to investigate any discordance between glycemic marker results.35 An age-stratified random sample of the population aged ≥18 years was drawn from Spain’s National Health System Registry. From November 2012 until March 2015, all subjects were successively convened for one day at the A Estrada Primary Care Centre for an evaluation which comprised fasting venous blood sampling, questionnaire interviews, and the description of subjects’ lifestyles. FPG was determined in subjects’ plasma samples using 5934 LADO-BALEATO et al. TABLE 1 Percentage coverage of the estimated conditional bivariate region for different sample sizes (500, 1000, and 2000), and nominal levels (5%,50%,90%and, 95%) in the three simulation scenarios X1 Sample size Nominal Global 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Scenario 1 n =500 5 4.6 4.4 4.8 4.7 4.7 4.8 4.7 4.6 4.6 4.7 4.7 4.3 50 47.9 46.7 48.5 48.3 48.3 48.6 48.6 48.4 48.4 48.2 48.0 45.3 90 88.3 87.3 88.8 88.6 88.7 88.8 88.8 88.7 88.7 88.4 88.2 86.0 95 93.6 92.9 93.9 93.9 94.0 93.9 94.0 93.9 93.9 93.7 93.5 91.9 n=1000 5 4.7 4.7 4.8 4.8 4.7 4.7 4.6 4.7 4.6 4.6 4.7 4.3 50 48.8 49.1 49.3 49.4 49.1 49.0 49.0 48.8 48.7 48.5 48.6 46.9 90 88.9 89.2 89.2 89.2 89.1 89.1 89.2 88.9 89.0 88.8 88.6 87.6 95 94.1 94.4 94.3 94.4 94.3 94.3 94.3 94.1 94.1 94.0 93.9 93.2 n=2000 5 4.7 4.7 4.8 4.8 4.7 4.8 4.7 4.7 4.7 4.7 4.8 4.6 50 49.2 48.9 49.7 49.4 49.3 49.4 49.5 49.3 49.3 49.1 49.1 48.0 90 89.2 89.1 89.3 89.5 89.3 89.4 89.5 89.3 89.2 89.4 89.1 88.4 95 94.3 94.3 94.3 94.6 94.5 94.4 94.5 94.4 94.4 94.4 94.2 93.7 Scenario 2 n =500 5 4.6 4.3 4.6 4.6 4.7 4.7 4.8 4.7 4.6 4.7 4.6 4.0 50 47.6 45.4 47.9 48.3 48.6 48.8 48.7 48.4 48.1 48.0 47.9 43.7 90 88.7 87.0 88.8 89.1 89.6 89.6 89.7 89.5 89.3 89.0 88.7 85.7 95 94.1 92.9 94.2 94.4 94.7 94.7 94.8 94.7 94.6 94.4 94.0 92.0 n=1000 5 4.5 4.4 4.6 4.5 4.6 4.6 4.5 4.6 4.5 4.5 4.4 4.3 50 47.0 46.1 46.9 47.1 47.3 47.7 47.4 47.3 47.5 47.3 47.0 45.9 90 88.5 87.3 88.1 88.4 88.7 88.9 88.8 89.0 88.9 88.8 88.8 88.3 95 94.0 93.0 93.6 93.8 94.1 94.2 94.2 94.4 94.3 94.3 94.2 93.9 n=2000 5 4.6 4.4 4.7 4.6 4.6 4.6 4.6 4.6 4.6 4.6 4.7 4.2 50 48.3 47.0 48.2 48.6 48.8 49.0 49.0 48.6 48.7 48.7 48.8 46.0 90 89.3 88.2 89.2 89.6 89.6 89.8 89.8 89.7 89.6 89.5 89.3 87.7 95 94.5 93.7 94.4 94.7 94.7 94.9 94.9 94.9 94.7 94.6 94.6 93.5 Scenario 3 n =500 5 4.6 4.4 4.7 4.7 4.8 4.7 4.7 4.7 4.6 4.6 4.5 3.9 50 47.0 44.9 47.2 47.7 48.2 48.4 48.2 47.9 47.7 47.5 46.9 42.3 90 88.3 86.1 87.9 88.7 89.2 89.4 89.4 89.3 89.1 88.7 88.3 85.6 95 93.8 92.1 93.4 94.0 94.4 94.5 94.5 94.5 94.4 94.1 93.8 92.1 n=1000 5 4.5 4.4 4.6 4.5 4.6 4.6 4.5 4.6 4.5 4.5 4.4 4.3 50 47.0 46.1 46.9 47.1 47.3 47.7 47.4 47.3 47.5 47.3 47.0 45.9 90 88.5 87.3 88.1 88.4 88.7 88.9 88.8 89.0 88.9 88.8 88.8 88.3 95 94.0 93.0 93.6 93.8 94.1 94.2 94.2 94.4 94.3 94.3 94.2 93.9 n=2000 5 4.9 4.4 4.6 4.5 4.6 4.6 4.5 4.6 4.5 4.5 4.4 4.3 50 49.5 46.1 46.9 47.1 47.3 47.7 47.4 47.3 47.5 47.3 47.0 45.9 90 89.6 87.3 88.1 88.4 88.6 88.9 88.8 89.0 88.9 88.8 88.8 88.3 95 94.6 93.0 93.6 93.8 94.1 94.2 94.2 94.2 94.3 94.3 94.2 93.9 LADO-BALEATO et al. 5941 where G(⋅)is a known link function, 𝛼is a constant, and fjunknown functions representing the effects of continuous covariates. Given a sample {(Yi,Xi)}n i=1, this model can be estimated using the following iterative process based on a Newton Raphson procedure which extends the ACE (Alternating Conditional Expectation) algorithm.50 Initialize: compute the initial estimates, 𝛼 =G−1(Y)with Y=n−1∑n i=1Yi, f0 1,…, f0 p=0. Step 1: for i=1,…,nconstruct the linearized response  Yand the weights Wso that:  Yi=𝜂0 i+Yi−G(𝜂0 i) G′(𝜂0 i)and Wi=G′(𝜂0 i)2 𝜎2 i , where 𝜂0 i=𝛼 +∑p j=1 f0 j(Xj),G′(𝜂)=𝛿G 𝛿𝜂 ,and𝜎2 iis an estimation of the variance 𝜎2(Yi|G(𝜂0 i)). The estimated 𝜎2 ican be obtained fitting an additive model to (Yi−G(𝜂0 i))2. Step 2: fit an additive model to  Yweighted by Wand compute the updates 𝛼 and  fjfor j=1,…,p.Atthisstepwe have used an inner backfitting algorithm based on a local polynomial kernel smoother: Step 2.0: update the constant 𝛼 =(∑n i=1Wi)−1∑n i=1Wi Y Step 2.1: for j=1,…,pcalculate the partial residuals Rj i= Yi−𝛼 − j−1 ∑ k=1  fk(Xik)− p ∑ k=j+1  f0 k(Xik) and for i=1,…,n, compute the polynomial kernel estimator updates:  fj(Xij)= 𝜓 (Xij,{(Xlj,Rj l,Wi)}n l=1,hf j) being hf jthe smoothing bandwidth associated with the estimation of fj. Step 2.2: Repeat Step 2.1 replacing  f0 j(Xij)by  fj(Xij)for j=1,…,pand i=1,…,n, until the convergence criterion ∑n i=1( fj(Xij)− f0 j(Xij))2 ∑n i=1( f0 j(Xij))2 +0.001 ≤𝜀for all j=1,…,p is reached. Step 3: repeat the Steps 1 and 2with 𝜂0 ibeing replaced by 𝜂i=𝛼 +∑p j=1 fj(Xij)for i=1,…,nuntil: |MSE(𝜂, Y)−MSE(𝜂0,Y)| MSE(𝜂0,Y)≤𝜖, where 𝜖is a small and the mean squared error MSE(𝜂,Y)is defined as MSE(𝜂, Y)=n−1∑n i=1Wi(Yi−G(𝜂i))2 The proposed algorithm use two loops: (i) an external loop, for adjusting the transformed response models (Steps 1 and 3), (ii) an internal loop which estimates the nonlinear effect of continuous covariates using a backfitting method (Step 2). This algorithm was used to estimate the additive models for means (Equation (5)) with identity link , the additive models for variances (Equation (6)) with exponential link and the correlation additive model (Equation (7)) with tanh(⋅)link. Note that in the first case the algorithm is reduced to the internal loop. Finally, the polynomial kernel smoother used in Step 2 may be replaced by an alternative estimator (eg, penalized splines).51 APPENDIX B. LOCAL POLYNOMIAL KERNEL SMOOTHERS Given a sample {(Xi,Yi)}n i=1with a vector of weights {Wi}n i=1the local linear kernel smoother at a location x,𝜓(x)= 𝜓 (x,{(Xi,Yi,Wi)}n l=1,h)is defined as 𝜓(x)=  𝜷,where  𝜷=( 𝛽0,…, 𝛽q)is a vector which minimizes: 5942 LADO-BALEATO et al. n ∑ i=1 Wi[Yi− q ∑ j=0 𝛽j(Xi−x)j]2 K(Xi−x h), where K(⋅)denotes a kernel function (a symmetric density), qis the polynomial degree, and h>0 is the smoothing parameter, chosen using: CV =1 n n ∑ i=1 Wi(Yi−𝜓(−i)(Xi))2, where 𝜓(−i)(Xi)indicates the fit at Xileaving out the ith data vector. APPENDIX C. BIVARIATE KERNEL BANDWIDTH SELECTION In this section, the bivariate kernel bandwidth estimator proposed in Equation (10) is evaluated and compared with the well-known plug-in, and least-square cross-validation estimators, and an arbitrary large bandwidth (5H). The influence of the kernel bandwidth on the shape of the final conditional reference region is determined by examining the RMSEs (see Equation (11)), and the region perimeter with respect to the theoretical regions. Moreover, the coverage properties of the estimated region is assessed using an out-sample design. This evaluation was performed under two scenarios: •Scenario 1: involving a bivariate response simulated from the modelstructure (Equation (1)). The bivariate error (𝜀1,𝜀 2) was drawn from a bivariate standard Gaussian distribution. •Scenario 2: involving a bivariate response simulated from a parametric copula representation (see Equation (13)), considering a reverse Gumbel and logistic distribution joined by a Gumbel copula. In both scenarios the bivariate response parameters were made dependent on a continuous (X1∈U[0,1]) and a binary regressor (X2), specifically: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝜇1(X)=1+X2+f1(X1) 𝜎2 1(X)=1+X2 𝜇2(X)=1+0.5X2+f2(X1) 𝜎2 2(X)=1+0.5X2 𝜌(X)=0.3+0.2X1+0.3X2 , where f1(X1)=X1sin(3X1)and f2(X1)=sin(2𝜋X1)represent nonlinear effects. The sample size were set to n=500,1000, and 2000 and the evaluation performed using 1000 replicates. Figure C1 shows the main results for both simulation scenarios. With Gaussian data (scenario 1), the RMSE of the estimated bivariate reference region (𝜏=0.95) is higher for the plug-in, and cross-validation methods for every sample size, and continuous covariate value. Moreover, estimated regions’ perimeters were higher than the real one for classical bandwidth estimators. Accordingly, a kernel bandwidth wider than plug-in would obtain a better estimate of the bivariate region. Given that the bivariate error is represented using a standard bivariate Gaussian distribution, which is also used as a kernel function in the density estimator, this result is not surprising. This also explains the good coverages obtained for a large kernel bandwidth (see Table C1). Hence, if the standardized bivariate residuals are Gaussian-distributed, a parametric expression would be better used in the estimation of R𝜏. For non-Gaussian data (scenario 2), our proposed method obtains a smaller estimation error than that returned by the plug-in, and cross-validation estimators. Moreover, the use of a large bandwidth returns a worse estimation error. In addition, theoretical region perimeters are better approximated by our kernel bandwidth estimator. Figure C2 shows the estimated bivariate reference region for every kernel bandwidth. Plug-in, and cross-validation estimators returned greater variabilities in their estimates, while the large bandwidth ignored the true shape of the region. Moreover, the best coverages were obtained using the bandwidth selector described in Section 2.3 (see Table C1). Specifically, the plug-in, and cross-validation methods seems to overfit the training dataset, resulting in a coverage below 95% for the estimated reference region, while the large bandwidth (5H) offer a coverage of LADO-BALEATO et al. 5943 n = 500 n = 1000 n = 2000 Perimeter RMSE 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 155 205 255 305 355 405 0.0 0.1 0.2 0.3 0.4 0.5 Covariate Value Best−Coverage Plug−in CV 5H Scenario 1: Gaussian data n = 500 n = 1000 n = 2000 Perimeter RMSE 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 240 290 340 390 440 490 540 590 640 690 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Covariate Value Best−Coverage Plug−in CV 5H Scenario 2: non−Gaussian data FIGURE C1 Bivariate reference region performance depending on kernel bandwidth estimator for different X1predictor variable values, with X2fixed at zero, sample sizes (500, 1000, and 2000), for Gaussian and non-Gaussian data. Red line represents the theoretical region’s perimeter. Best-Coverage represents kernel bandwidth estimator proposed in Equation (10), CV is the least-square cross-validation method, and 5Han arbitrary large bandwidth [Colour figure can be viewed at wileyonlinelibrary.com] over 95%. Hence, when dealing with nonstandard responses, the density estimator plays a key role in the performance of the bivariate reference region. The present method shows better performance than the plug-in, and cross-validation estimators, but the bandwidth selection problem requires further work. APPENDIX D. BOOTSTRAP INFERENCE In this section we present a bootstrap procedure to obtain punctual confidence intervals, given a specific vector of covariates X0, for the components (mean, deviation and correlation components ) of the model presented in Equation (1). The steps for construction of the bootstrap confidence intervals are: Step 1. From the sample data {(Yi1,Yi2),Xi}n i=1obtain the estimates 𝜇r(X0),𝜎r(X0)(r=1,2)and 𝜌(X0). 5944 LADO-BALEATO et al. Best Coverage Plug−in Cross−validation 5H FIGURE C2 Estimation of the bivariate reference region for 50 replicates (gray) along with the theoretical region (red), for n =1000, X1=0.5andX2=0, and the bivariate kernel bandwidths used in the estimation of the model’s residual density function. Best coverage represents our Equation (10) proposal, and 5Han arbitrary large bandwidth [Colour figure can be viewed at wileyonlinelibrary.com] Step 2.Forb=1,…,Bgenerate bootstrap samples {(Y• i1,Y• i2),Xi}n i=1with (Y• i1 Y• i2)=(𝜇1(Xi) 𝜇2(Xi))+ 𝜮1∕2(Xi)(𝜀• i1 𝜀• i2), where {(𝜀• i1,𝜀• i2)}n i=1is a sample of size nfrom the residuals {(𝜀i1,𝜀i2)}n i=1with replacement, and compute 𝜇•b r(X0),𝜎•b r(X0) and 𝜌•b(X0)as in Step 1. The limits for the 100(1−𝛼)% confidence intervals of the true components 𝜇r(X0),𝜎r(X0),and𝜌(X0)are given respectively by (𝜇𝛼∕2 r(X0),𝜇1−𝛼∕2 r(X0)),(𝜎𝛼∕2 r(X0),𝜎1−𝛼∕2 r(X0)),and(𝜌𝛼∕2(X0),𝜌1−𝛼∕2(X0)),where 𝜇p r(X0)represents the p-percentile of 𝜇•1 r(X0),…,𝜇•B r(X0),𝜎p r(X0)represents the p-percentile of 𝜎•1 r(X0),…,𝜎•B r(X0),and 𝜌p(X0)is the p-percentile of 𝜌•1(X0),…,𝜌•B(X0). APPENDIX E. R CODE AND DATA ANALYSIS The statistical methods developed in this article, and the AEGIS dataset, are available in the R package refreg.52 In the following we present the R code necessary to replicate Section 4.2 results. install.packages("refreg") library(refreg) #-- AEGIS dataset ?aegis head(aegis) summary(aegis) dm_no =subset(aegis, aegis$dm=="no") # healthy patients sample #-- Fitting bivariate location-scale model mu1 =fpg ∼s(age) LADO-BALEATO et al. 5945 TABLE C1 Coverage probability of bivariate data points for different sample sizes (500, 1000, and 2000), covariate X1 values, with X2fixed at zero, and bivariate kernel bandwidths estimators X1 Sample size Bandwidth 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Scenario 1 n =500 Best coverage 93.0 93.9 93.9 94.0 94.0 94.0 93.8 93.7 93.6 93.4 91.7 Plug-in 91.5 92.5 92.6 92.6 92.7 92.6 92.5 92.3 92.3 92.1 90.1 CV 87.2 88.4 88.5 88.5 88.6 88.5 88.4 88.2 88.1 87.9 85.6 5H93.6 94.4 94.4 94.5 94.5 94.5 94.3 94.2 94.2 94.0 92.3 n=1000 Best coverage 93.8 94.3 94.3 94.3 94.3 94.3 94.2 94.1 94.1 94.0 93.2 Plug-in 93.1 93.6 93.6 93.6 93.6 93.6 93.5 93.4 93.4 93.3 92.4 CV 89.8 90.5 90.5 90.5 90.5 90.5 90.4 90.3 90.3 90.2 89.0 5H94.3 94.7 94.7 94.7 94.7 94.7 94.6 94.6 94.5 94.4 93.7 n=2000 Best coverage 94.2 94.5 94.5 94.5 94.5 94.4 94.4 94.4 94.3 94.3 93.8 Plug-in 93.8 94.1 94.1 94.1 94.1 94.1 94.0 94.0 93.9 93.9 93.4 CV 91.6 91.9 91.9 91.9 91.9 91.9 91.8 91.8 91.7 91.7 91.1 5H94.5 94.8 94.8 94.8 94.8 94.8 94.7 94.7 94.7 94.6 94.2 Scenario 2 n =500 Best coverage 92.5 93.6 94.0 94.3 94.3 94.2 94.0 94.2 93.6 93.4 91.7 Plug-in 89.9 91.1 91.8 92.0 92.2 92.2 92.0 92.1 91.5 91.1 88.4 CV 86.1 87.5 88.3 88.6 88.8 89.0 88.7 88.7 88.2 87.6 84.2 5H94.0 94.8 95.2 95.4 95.3 95.2 94.9 95.1 94.5 94.4 93.5 n=1000 Best coverage 92.7 93.5 93.9 94.0 94.0 94.0 94.1 94.4 94.0 94.1 93.5 Plug-in 91.2 92.1 92.5 92.7 92.9 92.9 93.0 93.4 92.9 92.9 91.6 CV 87.4 88.5 89.1 89.3 89.5 89.6 89.7 89.9 89.6 89.5 87.6 5H94.2 95.0 95.2 95.2 95.1 95.0 95.0 95.4 95.0 95.1 95.3 n=2000 Best coverage 93.1 93.6 93.9 94.1 94.1 94.1 94.2 94.6 94.1 94.1 93.7 Plug-in 92.2 92.7 93.0 93.3 93.3 93.5 93.5 93.9 93.4 93.4 92.8 CV 90.0 90.5 91.0 91.2 91.3 91.5 91.5 91.9 91.4 91.4 90.6 5H94.4 94.9 95.1 95.2 95.1 95.1 95.0 95.4 94.9 95.0 95.1 Note: Best coverage represents our Equation (10) proposal, CV is the least-square cross-validation method, and 5Han arbitrary large bandwidth. mu2 =hba1c ∼s(age) var1 =∼s(age) var2 =∼s(age) rho =∼s(age) formula =list(mu1,mu2,var1,var2,rho) fit =bivRegr(formula,data=dm_no) s_b =summary_boot(fit,B=250,parallel =T) # bootstrap 95% CI plot(s_b,eq=1) plot(s_b,eq=2) plot(s_b,eq=3) plot(s_b,eq=4) plot(s_b,eq=5,select=1)+ggplot2::geom_hline(yintercept=0,color ="red") 5946 LADO-BALEATO et al. #-- Bivariate reference region region0 =bivRegion(fit,H_choice ="plug.in",tau=c(0.50,0.90,0.95)) plot(region0,tau=0.95,col="gray",reg.lwd=2) # overfitted region set.seed(719) region1 =bivRegion(fit,H_choice ="Hcov",tau=0.95,k=50) # best coverage region plot(region1,tau=0.95,col="gray",reg.lwd=2.75, legend=T,xlab="Fasting Plasma Glucose, mg/dL", ylab="Glycated hemoglobin, %" # standarized region) summary(region1) # identify outside patients #-- Depict conditional region plot(region1,cond=T,newdata=data.frame(age=c(20,30,40,50,60,70,80)),col=NA, reg.lwd=2)