Geographical differences in blood potassium detected using a structured additive distributional regression model

Geographical differences in blood potassium detected using a structured additive distributional regression model Jenifer Espasandín-Domínguez a,*, Alfonso Javier Benítez-Estévez b, Carmen Cadarso-Suárez a, Thomas Kneib c, Tegra Barreiro-Martínez b, Balbina Casas-Méndez d, Francisco Gude e a Group of Biostatistics, Center for Research in Molecular Medicine and Chronic Diseases (CiMUS). University of Santiago de Compostela, Spain. b Department of Laboratory Medicine. Hospital Clínico Universitario de Santiago de Compostela, Spain. c Chair of Statistics. Georg-August-University Göttingen, Germany. d Department of Statistics, Mathematical Analysis, and Optimization. University of Santiago de Compostela, Spain. e Clinical Epidemiology Unit. Hospital Clínico Universitario de Santiago de Compostela, Spain. * Correspondence to: Group of Biostatistics, Center for Research in Molecular Medicine and Chronic Diseases (CiMUS), Campus Vida, 15782, Santiago de Compostela, A Coruña, Spain. E-mail address: [email protected] (J. Espasandín-Domínguez). Keywords: Potassium; Distributional regression; Spatial análisis; P-splines. 2 1. Introduction Recently, general practitioners working in the Santiago de Compostela Health Area (SCHA) in northwestern Spain raised concerns over the high percentage of patients whose serum potassium concentrations were above the normal range, and over differences in the values recorded from one area to another. Analytical laboratories are commonly called upon to determine serum potassium concentrations, especially for patients with diabetes, heart and kidney disease. When potassium concentrations are recorded falsely as high (pseudohyperkalaemia) owing to specimen-collection or processing errors, medical mistakes can be made with disastrous consequences for patients. Although the list of sample management factors that can modify the potassium concentration is large, problems can be prevented by good laboratory practice (Stankovic and Smith, 2004). The SCHA covers an area of approximately 4905 km2, and at the time of the present study had a population of 497 171 (Instituto Galego de Estatística, 2015). The Clinical and Laboratory Standards Institute (2008) recommends that procedures be established for the transport of samples to laboratories to ensure that they are protected from deterioration (Tanner et al., 2008). The 46 general practices that use the laboratory service based at the Hospital Clínico Universitario de Santiago (CHUS) are, however, up to 70 km away, and timely transport of samples to the laboratory is a challenge. The aim of the present work was to determine whether any geographical differences exist in terms of recorded serum potassium concentrations and their variability that might be attributed to preanalytical factors, such as the centre where blood was extracted, adjusting for other potential covariables that might influence the results. For this, a structured additive distributional regression model (Klein et al., 2015) was used. This type of model is closely related to generalized additive models for location, scale and shape (GAMLSS, Rigby and Stasinopoulos, 2005). Inference in GAMLSS models is based on penalized maximum likelihood estimation, achieved via backfitting loops over the additive predictor components. The current work present a Bayesian version of a structured additive distributional regression model dependent on Markov chain Monte Carlo (MCMC) simulation algorithms. This approach has the advantage of providing credible intervals without relying on asymptotic arguments (Klein et al., 2015). An advantage of this kind of model is the possibility of incorporating spatial effects. However, in most cases the output of spatial effects is not directly interpretable by biomedical researchers. This paper proposes a way in which spatial effects can be visualized. Another advantage of this type of model is the possibility of contemplating a wide range of response variables. The deviance information criterion (DIC) (Spiegelhalter et al., 2002;Klein et al., 2015) is commonly used for model choice in distributional regression. Quantile residuals can be used to check the performance of a selected model (Klein et al., 2015). In practice, the residuals can be assessed graphically in terms of quantile–quantile plots. However, interpreting the resulting graphs can be difficult, and the decision on the adequacy of a model remains subjective. Sometimes, even though the model is correct, the plot may deviate substantially from a straight line (Augustin et al., 2012). We therefore here propose the use of quantile–quantile plots with reference bands. To construct these bands, the methodology of Augustin et al. (2012) was adapted to the context of distributional regression. In Section 2we present the description of the database used in the study. Section 3provides an introduction to structured additive distributional regression models, including the choice of the response distribution (and thus the model selected) and the construction of quantile–quantile plots ABSTRACT: Recently, physicians in an area of northwestern Spain became concerned about the large number of patients whose serum potassium concentrations were above the normal range, as well as differences in the values recorded from one area to another. With the aim of identifying geographical differences in both mean and variability of potassium levels, analyses were performed using modern flexible regression techniques based on a structured additive distributional regression model. In this type of model, every parameter of a response distribution – rather than just the mean – is related to a structured additive predictor. After adjusting for variables such as age, sex, clot-contact time and spatial effects, differences in potassium concentrations were confirmed. The type of distributional regression model used permitted the mean and variance of the potassium concentrations to be modelled using additive predictors that allow for different types of covariate effects. A variety of complex distributions were contemplated. In general, higher concentrations and greater variability were recorded in areas further from the hospital laboratory responsible for the analyses, although some that were nearby also returned high values. Further actions are required to confirm whether these differences reflect reality or non-optimum pre-analytical handling (resulting in pseudohyperkalaemia) in certain districts. 3 with reference bands. Section 4provides the results obtained following analysis of the potassium database discussed in Section 2. Finally, Section 5provides a summary and some comments on directions of future research. 2. Data description The database used in this work was provided by the Clinical Analysis Laboratory of the CHUS. This supplied information on all blood extractions performed between 1 June and 31 December 2015 for which serum potassium, sodium and creatinine measurements were made. The initial number of samples was 145 960, collected at 46 extraction centres within the SCHA. Those samples showing signs of haemolysis were excluded, as were those with creatinine or sodium concentrations outside the normal range (indicators of impaired kidney function). The final number of samples used in the present analysis was therefore 95 096. 2.1. Health Area of Santiago de Compostela (SCHA) The SCHA, in Spain’s northwest, covers 46 municipal districts (see Table 2 and Fig. 6 in the Appendix for a map and the distribution of the population). The CHUS reference hospital is located in the city of Santiago de Compostela, from where the SCHA’s health centres and doctors’ practices are coordinated (see Fig. 6 in the Appendix for their locations). Blood samples were taken, usually daily, at designated locations in these 46 areas and transported by road to the CHUS following different routes. The following variables were considered to be covariates: gender,age (in years), clot-contact time in minutes (cctime in formulae (1)), and demographic information on the district where the extraction centres were located (s). 56% of the patients who provided blood samples were female, and 44% were male. The age range of the patients was 1–103 years. The mean (respectively SD) age was 54.8 (19.0) years; the median was 61 years. The clot-contact time was taken as the difference between the starting time of sample collection at the extraction centre and the entry time of the sample in the laboratory registry. The range was 4–458 min, the mean (SD) clot-contact time 229 (49) min, and the median time 231 min. 3. Structured additive distributional regression models Distributional regression models (Klein et al., 2014) are innovative models that permit marginal distribution parameters to be modelled using additive predictors that allow for several types of covariate effects (such as the non-linear effects of continuous covariates, random effects, and the interactions or spatial effects). By modelling each parameter of the response at the same time – and not just the mean – they provide additional flexibility. This type of model also provides different types of (possibly non-standard) response distributions for continuous, discrete, and mixed discrete continuous distributions. However, the adequate selection of response variable in the formulation of these models is not without its difficulty. The statistical literature usually suggests one suppose the response variable to follow a normal distribution. However, in the present work, several distribution types were contemplated for the modelling of the potassium concentrations, including the log-normal distribution, the truncated normal distribution, the inverse Gaussian distribution, and the gamma distribution (see Section 3.3 and Appendix A.2). In this work, a structured additive distributional regression model was used to study the potassium concentrations recorded. This type of distributional regression allows the effect of the covariate information on all the parameters of the response distribution to be examined. Let us assume that observations (yi,νi,i=1, . . . , n) are made, where yiare observations on the response variable, and νirepresents the generic covariate vector. In this scenario, the response variables yican be assumed independently distributed with K-parametric densities p(yi|ϑi1, . . . , ϑiK )≡pi. In other words, the conditional distribution piof an observation yigiven νiis expressed in terms of the Kdistributional parameters of the response distribution: ϑi1, . . . , ϑiK . 4 In structured additive distributional regression models, each parameter ϑik,k=1, . . . , K, of the response distribution is related to a semiparametric additive predictor ηϑk idefined in terms of the covariate vector νi. As in other types of classic regression model, such as generalized linear regression models, a suitable response function is used to map the predictor to the parameter of interest, ϑik =hϑk(ηϑk i). Following Klein et al. (2015), in this expression, the superscript ϑkrefers to the fact that Kpredictors specific, for each of the distribution parameters of the response variable (and not just for the mean as in classical regression), are taken into account. Moreover, for an observation i=1, . . . , n, a suitable structured additive predictor for parameter ϑkcan be written as: ηϑk i=βϑk 0+fϑk 1(νi)+···+fϑk Jk(νi)+fϑk spat (si) where βϑk 0represents the overall level of the predictor, and the functions fϑk j(νi),j=1, . . . , Jk represent the different covariate effects. Note that each distribution parameter may depend on different covariates and a different number of effects, say Jk. The generic representation with the complete covariate vector can be used to simplify the notation. Finally, fspat (s) is the spatial effect capturing heterogeneity at the level of the districts s. In structured additive regression, each function fjis approximated by a linear combination of Dj appropriate basis functions: fj(νi)= Dj ∑ dj=1 βj,djBj,dj(νi). In matrix notation, we can write fj=(fj(ν1), . . . , fj(νn))′=Zjβjwhere Zj[i,dj] = Bj,dj(νi) is an n×Djdesign matrix and βjis the vector of coefficients (with dimension Dj) to be estimated. The basis function representation then lead us to the following matrix representation of the predictor 3: η=β01+Z1β1+···+ZJβJ. For each of the parameter vectors βj, the multivariate normal prior can be assumed: p(βj|τ2 j)∝(1 τ2 j) rk(Kj) 2 exp (−1 2τ2 j β′ jKjβj) in which the (potentially rank-deficient) precision matrix Kjcorresponds to the penalty matrix in a frequentist formulation. Note that here, we are using a generic notation for different terms (e.g penalized splines, Markov random fields, random effects). The precision matrix, Kj, for the Markov random field is rank-deficient by construction since only deviations from a constant spatial effect are penalized. This leads to a rank deficiency of one. In the case of the penalized splines, the rank deficiency comes from the fact that polynomials of degree equal to the difference order minus one are not penalized. However for other possible terms like random effects this matrix is not rank-deficient. A prior smoothing variance of τ2 jis assigned as an inverse gamma hyperprior τ2 j∼IG(aj,bj) (with aj=bj=0.001 as a default option) in order to obtain data-driven smoothness. Again following Klein et al. (2015) and to simplify the notation, the dependence on the distributional parameter indicated by the superscript ϑk, the observation index i, and the function index j, have been dropped. Fahrmeir et al. (2013) discuss all terms included in this generic predictor. The following paragraphs outline the prior assumptions for the hierarchical predictor required for this study. See Lang and Brezger (2004) for more details. 3.1. Linear effects and continuous covariates For the effect of the intercept, β0, and the gender of the individuals, β1, a flat, non-informative prior was assumed. 5 The non-linear effects of continuous covariates (age and clot-contact time) were modelled using Bayesian versions of penalized splines (P-splines, Lang and Brezger, 2004), introduced into a frequentist setting by Eilers and Marx (1996). To model age and clot-contact time, 20 inner knots, a cubic spline basis, and a second order random walk prior for penalized splines were contemplated. For the penalized spline specifications, we were able to rely on extensive research concerning the number and placement of knots, the order of the random walk prior, and the degree of the polynomial spline, e.g. Eilers and Marx (1996), Lang and Brezger (2004) and Brezger and Lang (2006). Their main findings can be summarized as follows: (i) The number and placement of the knots has only a very minor impact on the fit if the number of knots chosen is not too small. (ii) 20 equidistant knots yield sufficient flexibility for basically all situations of applied interest. (iii) Second order random walk priors leave a linear effect unpenalized which is in analogy to the common penalty for smoothing splines. Moreover, first order differences often yield more wiggly estimates. (iv) Finally, cubic splines yield a visually smooth function estimate which is twice continuously differentiable. This fits very well with the common visual perception of non-linear effects. 3.2. Spatial effects The spatial effects, fspat , were understood as the sum of the spatially structured correlated (smooth) effects fstr , and spatially uncorrelated (unsmooth) effects: fspat (s)=fstr (s)+funstr (s). A spatial effect is usually a surrogate of many unobserved influential factors, some of which may obey a strong spatial structure while others may be present only locally. By estimating a structured and an unstructured component, it was hoped that distinctions could be made between these kinds of influential factor (Besag et al., 1991). For correlated spatial effects (or structured spatial effects), we assume spatial correlations defined implicitly by assuming a Markov random field (Fahrmeir et al., 2013) as a prior distribution for the separate regression coefficients corresponding to the distinct regions. The Markovian structure is determined by the neighbourhood structure for the regions and the precise form of the prior distribution is defined by: βstr,s|βstr,r,r= s, τ2 str ∼N(1 Ns∑ r∈δs βstr,r,τ2 str Ns), where Ns=|δs|is the number of adjacent sites or neighbours, and r∈δsdenotes that region ris a neighbour of site s. The conditional mean of βstr,s, given all other coefficients, is the average of the neighbouring regions. Additional uncorrelated random effects (or unstructured spatial effects) may be incorporated as a surrogate for unobserved local small-area, group or individual specific heterogeneity. If one ignores spatial proximity and interprets s∈{1, . . . , S}as a cluster variable that only represents membership to different groups (such as different individuals in longitudinal data or more generally to different clusters of observations), we can assume a standard Gaussian random effects prior, i.e. funstr (s)∼ N(0, τ2 unstr ),s∈{1, . . . , S},where the different groups correspond to the different administrative regions in the data set. Note that here two variances are used for the structured, τ2 str , and the unstructured effect, τ2 unstr . Structurally, both variances are of the same type but they refer to different prior assumptions, where one assumes spatial structure (τ2 str ) while the other one assumes spatial independence (τ2 unstr ). 3.3. Inference and choice of the response distribution Distributional regression models can be inferred employing computationally efficient extensions of the MCMC techniques developed by Klein et al. (2014). Many approaches have been proposed for dealing with model choice in a Bayesian framework. In the present work, the Deviance Information Criterion (DIC, Spiegelhalter et al., 2002;Klein et al., 2015) 6 Fig. 1. Quantile–quantile residuals plot for the selected model with reference bands: the closer the residuals to the bisecting red line, the better the fit to the data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) was used to choose the best response distribution. The DIC is a commonly used criterion for model choice in Bayesian inference. It became popular in part because of its easy implementation from the MCMC output. The performance of the DIC was valued as positive by Klein et al. (2015), who compared several mis-specified models with the true DIC model. In the present work, the DIC showed the lognormal distribution to provide the best and most parsimonious fit (see Appendix A.2 for more details). The residuals can also be used to check the performance of a selected model (Klein et al., 2015). If the estimated model is close to the true model, the quantile residuals approximately follow a standard normal distribution, even if the model distribution itself is not a normal distribution. In practice, the residuals can be assessed graphically in terms of quantile–quantile plots. To improve the interpretability of this type of graphics, we propose plotting reference bands around the diagonal line to give a rough indication of the uncertainty implied by estimating the model from finite data. More precisely, following Augustin et al. (2012), we repeatedly simulate data from the fitted model and add pointwise minima and maxima for the resulting quantile residuals from a pre-specified number of replications. Fig. 1 shows quantile–quantile plots for the selected model with a log-normal distribution. Fig. 1 also provides reference bands for judging the relevance of departures of quantile–quantile plots from the ideal red line. The log-normal distribution turns out to be appropriate for residuals in the range between −2.5 and 2.7 but deviates from the diagonal line for extreme values. Note, however, that these extreme values correspond to only 1.37% of the total data (0.87% of the database to the left and 0.5% to the right) such that our model explains the vast majority of observations well. Thus, the structured additive distributional regression model described at the beginning of Section 3was used with a log-normal distribution response and with two covariate-dependent parameters (corresponding to the mean of the log-transformed potassium concentrations and the scale parameter σ2), to describe and compare the potassium concentrations recorded in the different districts of the SCHA. This model is expressed as follows: {ηµ=βµ 0+gender′βµ 1+fµ 1(age)+fµ 2(cctime)+fµ spat (s), ησ2=βσ2 0+gender′βσ2 1+fσ2 1(age)+fσ2 2(cctime)+fσ2 spat (s) (1) where β0represents the overall level of the predictor and β1captures the effect of the gender. Moreover, f1(age) and f2(cctime) are non linear effects of the age and the clot-contact time, respectively. Finally, fspat (s) is the spatial effect of the districts s. 4. Data analysis Statistical analyses were performed using open-source BayesX software (Belitz et al., 2015). The BayesX (Kneib et al., 2015) and R2Bayesx (Umlauf et al., 2015;Belitz et al., 2016) R packages were used as graphic interfaces. 7 Table 1 Summary of estimated linear effects for model (1). Parameter mean 2.5% quantile median 97.5% quantile βµ 0(intercept) 1.489 1.465 1.484 1.512 βσ2 0(intercept) −0.032 −0.050 −0.031 −0.014 βµ 1(gender) −0.006 −0.008 −0.006 −0.005 βσ2 1(gender) −0.032 −0.050 −0.032 −0.014 (a) fµ 1(age). (b) fσ2 1(age). Fig. 2. Posterior mean estimates of non linear effects of age on µand σ2. (a) fµ 2(cctime). (b) fσ2 2(cctime). Fig. 3. Posterior mean estimates of non linear effects of cctime on µand σ2. All results are summarized in Table 1 and Figs. 2–5. For the continuous variables, all covariates but the one that is visualized are fixed at their average while spatial effects are set to zero. For the spatial effects, we used district-specific averages for all covariates and determined significance based on the comparison with the average of all spatial effects. 8 Fig. 4. Posterior mean estimates of the complete spatial effects on mean potassium levels, fµ spat ,and 95% posterior probabilities (right). In the right panel a value of 1 corresponds to a strictly positive 95% credible interval, and a value of −1 to a strictly negative credible interval; a value of 0 indicates that 0 is contained in the corresponding credible interval. Fig. 5. Posterior mean estimates of the complete spatial effects on the variance of potassium levels, fσ2 spat , and 95% posterior probabilities (right). In the right panel a value of 1 corresponds to a strictly positive 95% credible interval, and a value of −1 to a strictly negative credible interval; a value of 0 indicates that 0 is contained in the corresponding credible interval. Gender had some influence on the results (both on the mean and variance of the potassium concentrations), but this effect is clinically not relevant since the differences between both men and women were minimal (0.04 mg/dL). The same for the variance. However, age (Fig. 2), clot-contact time (Fig. 3) and the place of origin of the samples did influence the potassium concentrations recorded (Figs. 4 and 5). Children and the elderly had higher mean potassium concentrations than did patients of intermediate age; the values recorded for the elderly also showed greater variability (Fig. 2). These age-related findings might, however, be expected since venipunture is harder to perform in both children and the elderly, and both age groups show greater capillary fragility. This can lead to situations in which haemolysis occurs, releasing potassium from the red blood cells and increasing the recorded concentration for both age groups, as well as the variability of values recorded for the elderly. The potassium concentrations recorded were clearly not uniformly distributed over the study area. In general, higher potassium concentrations were recorded in the areas farthest away from the test laboratory, although some areas close to it also returned high values. In some areas, these high concentrations were accompanied by greater variability in the results, particularly in the districts on the northern periphery of the study area. The districts to the southeast also returned high potassium concentrations but with the less variability. Although caution should be exercised when interpreting these spatial analysis results, it may be that potassium concentrations in the periphery are related to pre-analytical factors associated with the extraction centres. The affected areas are also those with the lowest population densities; they are therefore likely to have less equipment, fewer personnel, 9 and perhaps less well trained personnel than in the more central districts. These periphery districts may also have older inhabitants, clinical practices may be less homogeneous, and they are the worst communicated with the test laboratory (Figs. 4 and 5). 5. Discussion Distributional regression models extends the use of generalized additive models (GAM, Hastie and Tibshirani, 1990) to situations in which the response distributions are non-standard, and in which not only the mean but multiple parameters are related to additive predictors via suitable link functions. Further, they allow additional flexibility by specifying structured additive predictors for each parameter of interest, and thus adjust for flexible non-linear effects of continuous covariates for which the smoothness is determined based on the data. They also allow the contemplation of spatial effects to capture unobserved spatial heterogeneity and spatial correlations, interaction terms such as varying coefficients or interaction surfaces, and cluster-specific random effects (Fahrmeir et al., 2013; Brezger and Lang, 2006). The use of a new structured additive distributional regression model allowed for the flexible modelling of the distribution of potassium concentrations with a potentially non-standard response type, and permitted covariate effects to be taken into account, including the smooth estimation of the effect of continuous variables, categorical covariates, random effects and possible spatial trends. Its use clearly identified differences in serum potassium concentrations among extraction sites after adjusting for other potentially influential factors, such as age, gender and clot-contact time. The spatial analysis revealed some districts to return higher mean serum potassium concentrations, and to show greater variability in terms of these results. Two geographically-related clusters were detected: (1) districts on the periphery of the study area that returned higher potassium concentrations (and showed greater variability in the results) than those in the central area, and (2) a number of districts that returned higher potassium concentrations independent of their location. Although classic regression analyses allow for the easy interpretation of results, they only focus on means, and may lead to erroneous conclusions when modelling complex data structures. The distributional regression models used in this manuscript provide a generic framework for performing regression analyses in which several parameters of a potentially non-standard response distribution are related to flexible regression predictors. This work shows how to visualize from a statistical viewpoint the results of distributional regression models in an analysis comprising spatial information. It is not sufficient to show the estimated spatial effects directly; rather the spatial effect has to be adjusted with respect to the covariates observed in the particular regions. We therefore plugged in covariate values obtained as spatially stratified averages for all other covariates and then compared against the overall mean of the spatial effect to determine significances. Another benefit of this type of model is the possibility of being able to contemplate a wide family of response variables. One way of examining the goodness of fit of the selected model is via quantile– quantile plots. However, conclusions drawn from such plots can be subjective. This manuscript proposes the use of quantile–quantile plots with reference bands. The present work examined different (non-standard) distributions that depended on the mean and variance of responses. Using the DIC, distributional regression with a log-normal response was used. This provided not only improved goodness of fit over classic distributions (e.g., a Gaussian distribution), but led to different results being obtained. Thus, although the non-linear effects of the covariables on the expectation were rather similar with both distributions, the log-normal distribution allowed differences to be identified in the variability of the spatial effects associated with the potassium concentrations that were undetectable when a Gaussian distribution was contemplated (data not shown). The majority of the central municipalities had larger populations, more health care personnel and more equipment, and followed protocols more strictly, which might explain the lower potassium concentrations they recorded (with both distributions) and their smaller variability (lognormal distribution). In medical organizations, examining clinical variation in medical practice is an important step to measuring efficiency and effectiveness in care delivery. From a statistical viewpoint, another important feature of this type of model, is the possibility of modelling the effects of the continuous covariables and spatial effects in a flexible, unified manner