Generalized Mixed Modeling in Massive Electronic Health Record Databases: what is a healthy serum potassium?

Converting electronic health record (EHR) entries to useful clinical inferences requires one to address computational challenges due to the large number of repeated observations in individual patients. Unfortunately, the libraries of major statistical environments which implement Generalized Linear Mixed Models for such analyses have been shown to scale poorly in big datasets. The major computational bottleneck concerns the numerical evaluation of multivariable integrals, which even for the simplest EHR analyses may involve hundreds of thousands or millions of dimensions (one for each patient). The Laplace Approximation (LA) plays a major role in the development of the theory of GLMMs and it can approximate integrals in high dimensions with acceptable accuracy. We thus examined the scalability of Laplace based calculations for GLMMs. To do so we coded GLMMs in the R package TMB. TMB numerically optimizes complex likelihood expressions in a parallelizable manner by combining the LA with algorithmic differentiation (AD). We report on the feasibility of this approach to support clinical inferences in the HyperKalemia Benchmark Problem (HKBP). In the HKBP we associate potassium levels and their trajectories over time with survival in all patients in the Cerner Health Facts EHR database. Analyzing the HKBP requires the evaluation of an integral in over 10 million dimensions. The scale of this problem puts far beyond the reach of methodologies currently available. The major clinical inferences in this problem is the establishment of a population response curve that relates the potassium level with mortality, and an estimate of the variability of individual risk in the population. Based on our experience on the HKBP we conclude that the combination of the LA and AD offers a computationally efficient approach for the analysis of big repeated measures data with GLMMs.


Background and motivation
Electronic Health Records (EHR) have been adopted in more than 80% of all medical practices in the United States.The contain information about vital signs, lab results, admissions, healthcare facility features (e.g.type, location and practice) and individual subject outcomes (e.g.death or hospitalizations).Due to the large number of patients, and the equally large number of features collected in each subject, HER "big data" are increasingly being mined in the hope that such analyses will generate useful clinical insights that enhance patient safety, improve health care quality and even manage costs [2,6,8,20].A unique feature of EHR data that complicates analytics is the repeated measures and unbalanced nature of the data, featuring multiple and unequal number of observations in individual patients (IP), seen by different healthcare practitioners (HCPs) across healthcare facilities (HCFs).However existing approaches to big data such as deep learning [17] completely ignore this feature of EHR data, missing on a unique opportunity to mine variation at the level of IP/HCPs or HCFs.Subject or healthcare facility focused inference requires the deployment of conventional Generalized Linear Mixed Models (GLMM) for the analyses of repeated measures at the group (IP/HCP/HCF) level.GLMMs can tackle the entire spectrum of questions that a clinical researcher would want to examine using EHR data.Such questions typically involve the analyses of continuous (e.g.biomarker values, vital signs), discrete (e.g.development of specific diagnoses), time to event (e.g.survival) or joint biomarker-discrete/time-event data.Notwithstanding these aspirations, the implementation of GLMMs in existing statistical environments (e.g.SAS, R) has been shown to scale poorly if it can execute at all execute at all in big datasets with more than a few thousand individuals [5,16].Application of GLMMs necessitates the evaluation of high-dimensional integrals, which even for the simplest EHR analyses involves tens of millions of dimensions (e.g. one for each patient).However, due to the "curse of dimensionality" numerical integration is no longer tractable either analytically or numerically in such high dimensional space.This is major roadblock that limits the applications of formal statistical methods to big datasets.Approximate methods for fitting large GLMMs within the computational constraints of standard multicore workstations [5,9,16] or even parallel architectures (e.g.LinkedIn's GLMix [23]) have been proposed in the literature.These methods tend to rely on specific features of the datasets being analyzed to speed up calculations e.g. through data partitioning, by applying the method of moments for estimation, parameter block updates or utilizing stochastic gradient descent.Empirical analyses [16] demonstrates that such approaches invariably trade statistical efficiency for speed, effectively discarding valuable information hidden in big data.Last but certainly not least, quantification of uncertainty about model estimates becomes extremely challenging theoretically and, in some cases, it is not addressed at all in publications (e.g.GLMix).The ideal approach to the analyses of EHR big data with GLMMs, would seek to eliminate or reduce the need for statistically inefficient approximations, while retaining the rigor, numerical precision and uncertainty quantification measures (e.g. standard errors) that one has come to expect and trust from analyses of small to medium sized data.In this work we show how these ideals can be maintained in big data analytics by implementing theoretically sound estimation methods for GLMMs in a computationally efficient manner.The theoretical background for achieving these goals is based on the Hierarchical likelihood (h-lik) approach to GLMM estimation [10][11][12].The h-lik was initially introduced as a theoretical framework for understanding the statistical properties of GLMMs but later received attention for numerical work.In the context of h-lik the Laplace approximation is used to replace multivariate integration by numerical optimization, paving the road for computationally efficient, potentially parallelizable implementations that are uniquely suited to tackle the challenges posed by EHR big data.However this potential is not realized in current implementations of the h-lik because of the reliance on the original, iterative estimating equation formulation which can be numerically challenging to implement [4].In this work, we circumvent these limitations by implementing h-lik method in the R package TMB [7].TMB combines the Laplace Approximation (LA) for the numerical evaluation of multivariate integrals, with Algorithmic Differentiation (AD) for the fast calculation of objective (log-likelihood) functions.TMB itself implements an approximate estimation method for non-linear statistical models (including GLMMs) [21,22] but to date has not been adopted for the much more theoretically rigorous, yet numerically more frail h-lik calculations.TMB provides the key missing ingredients for fast, scalable h-lik calculations: AD, parallel computations to speed up the most computationally challenging steps and built-in uncertainty quantification features.4) We quantify the interindividual and healthcare facility variation around the curve that relates the K + level to survival.Assessment of variability is a major advantage of GLMMs over other analytical approaches, e.g.deep learning and has implications for both policy and clinical practice.This paper is organized as follows.Section 2 reviews estimation techniques for GLMMs and the h-lik approach.We also describe our implementation of the hlik calculations for GLMMs using the TMB package in the language R. Section 3 introduces Health Facts and presents our analysis of risk in relation to the potassium level and renal function.In the same section we compare the results obtained by TMB and the referent GLMM implementation in R using Adaptive Gaussian Hermite (AGH).Since AGH does not scale favorably with dataset size, we use small random samples (1% and 10%) of the final dataset.In section 4 we discuss the subject matter significance of the analyses presented herein and we discuss the impact of the proposed computational methodology for EHR big data.The proposed computational methodology has applications outside the narrow field of EHR analytics and we discuss those, as well as the path towards parallel implementations last.

Methods and Technical Solutions 2.1 GLMMs and the H-likelihood :
Let   , i = 1 …  denote data and let   and   stand for vectors of features ("covariates") for each data element in the model.The   feature vector is multiplied by the global ("fixed effects") p dimensional coefficient vector  that models the overall response, while the   by the group specific ("random effects") q dimensional vector .In the simplest case, the random effects feature vector indicates the group membership of the data, e.g. which observations came from the same individual, but more complex scenarios are possible dependent on the generation mechanism and the resulting repeated measures structure of the data.A GLMM is defined by the following properties: 1.The probability density function of the data ( | ), conditional on the random effects and the two feature vectors, is a member of the exponential family.The nature of the data determines the exponential family, e.g. for binary classification this would be the binomial distribution, for the analysis of time-to-event or counts of events (as in the present case) the Poisson distribution [1], for continuous data the Gaussian distribution.2. The expectation of the data μ  conditional on the random effects is given by the equation: g(μ  ) =     +   T  .The function g(⋅) is the link function and is determined by the nature of the regression e.g. it is the logistic map for binary classification, or the logarithmic function when modeling count data.The variance of the data is determined by the conditional mean and the link function, e.g. it is equal to the mean in the Poisson case which concerns us here.3. The random effects vector is assumed to have a zero mean multivariate normal distribution and covariance matrix G = G(γ} where γ are the variance components, i.e.   ~N(0, ) Learning in a GLMM context corresponds to estimating the values of the variance components (γ), the fixed () and the random effects coefficient vectors ().The h-lik is the statistical procedure for effecting such learning.In the h-lik the starting point is the joint likelihood of the data i.e. the product of  | and   .The h-likelihood function is the logarithm of this product and is used for parameter estimation in a staged, hierarchical fashion: 1. Learning about the variance components is based on maximization of the Restricted Extended Likelihood (REML) criterion that integrates the joint likelihood over the "nuisance" random and fixed effect parameters: 2. Learning about β and  is based on the marginal likelihood obtained by conditioning ("pluggingin") of the estimate of the variance components  ̂ obtained in the previous step: Computationally the h-lik uses the LA during optimization.For a general p dimensional integral, the LA assumes the following form: where  ̂ is the optimum of the function ℎ(, ) and H(ℎ,  ̂) is the Hessian matrix of the second partial derivatives of the function ℎ(, ) with respect to u evaluated at the optimum.All implementation of the hlik approach to date [4,14,18] follow the original algorithm by Lee and Nelder [10,12] in which a nested optimization procedure is used to maximize over the variance components (outer) and the fixed effects (inner optimization).Typically, a quasi-Newton algorithm such as L-BFGS is used for these optimizations.Estimates of the random effects are obtained by the iterative solution of a large linear system.Due to this complexity, implementations of the h-lik are error prone and sensitive to technical decisions about matrix representation and the finite difference scheme used to calculate gradients during optimization.

H-likelihood via the Laplace Approximation and Algorithmic Differentiation
For a theoretically sound statistical learning process it is only required that inferences about the variance components be based on (1) and for the fixed and random effects on (2).The curse of dimensionality implied the multivariate integration is substantially mitigated by the LA used in h-lik, so that the major computational bottleneck center around the optimization of (3).This is precisely the type of problems addressed by the library TMB in the R language [7].To use this library, one specifies the joint log-likelihood as a C++ template function, which is converted by the package into a computational graph ("tape").During tape generation, explicit dependencies between the parameters are identified by TMB and C++ code is automatically generated for the evaluation of the gradient of the log-likelihood and its Hessian.AD thus comes with two key advantages: a) the gradients are computed with accuracy that rivals that of hand-coded analytical expressions b) optimized software code for the Hessians also becomes available to the analyst.Algorithmic gradients overcome both the slowness and the numerical instability of finite difference schemes that have been used in numerical optimizations by h-lik implementations.Having code for the Hessian allows for the use of Newton rather than Newton optimization and is also required for uncertainty quantification.For GLMM applications, the Hessian matrix will be a sparse one (i.e.only a few entries will be non-zero) and TMB uses a coloring algorithm to detect the sparsity pattern in each specific dataset analyzed.This effectively reduces the computational and storage complexity for the Hessian calculations.

Parallel computations:
In large datasets, the major computational bottlenecks involve the computation of the determinant of the Hessian e.g. through a sparse Cholesky decomposition in (3) and the evaluation of the joint log-likelihood.In our analyses we used the parallel implementation of the BLAS libraries (MKL® by Intel that are distributed with Microsoft R Open) as suggested by the creators of the TMB package.The evaluation of the joint log-likelihood was also done in parallel using OpenMP which is natively supported by TMB.

Implementing the h-lik in TMB:
Implementation involves steps in two languages: C++ and R. Using C++ macros, the user specifies the joint log-likelihood, receives parameters, random effects, data and initial values from R. R is used to prepare data, generate initial values, generate the tape, invoke the optimizer and post process the results.In particular, an analyst: 1. Specifies the joint log-likelihood in C++ and compiles it into a dynamic library.This is done only once, and the dynamic library may be used with different data for the same general model.Hence TMB code for GLMMs is reusable across application domains.2. Generates the computational tape for a specific dataset designating parameters as either fixed or random.Gradient code becomes available for use by the user at this step.3. Calls an optimizer of their choice (we typically use the BFGS or its limited memory variant) which repeatedly evaluates the joint-likelihood and its gradient until convergence.4. Runs the report function of TMB from within R to generate standard error of all estimates in order to quantify uncertainty upon convergence.
When implementing the h-lik steps 2 and 3 must be called twice.First, we integrate over ,  by designating both as "random" when generating the tape.After optimization we obtain a point estimate for the variance components.Then fixing the value of the variance components, one generates the tape a second time designating  as random.After optimization, one calls a report to obtain parameter and uncertainty estimates for the fixed and random effects.

3.1
Health Facts: Cerner Health Facts is a comprehensive source of de-identified, real-world data that is collected as a by-product of patient care from over 700 healthcare facilities across the United States.The relational database contains clinical records with timestamped information on pharmacy, laboratory, admission and billing data for over 69 million unique patients.Types of data available include demographics, encounters, diagnoses, procedures, lab results, medication orders, medication administration, vital signs, microbiology, surgical cases, other clinical observations, and health systems attributes.Detailed pharmacy, laboratory, billing and registration data go as far back as 2000 and contain more than 630 million orders for nearly 3,500 drugs by generic name and brand.The two largest tables in the database contain more than 4.3 billion laboratory results and 5.6 billion clinical events linked to more than 460 million patient visits.For the purpose of our analyses we used the entire Health Facts record abstracting potassium measurements that had been obtained within 24 hours of a clinical encounter.Two subsets involving 1% and 10% of patients who had at least two potassium measurements were also analyzed.Since 50% of all patients had one only one measurement, these datasets were only 0.5% and 5% of the full dataset.Sub-setting the individuals with two measurements was done because AGH quadraturebased implementations which were used in comparisons (see below) often failed to converge without more than one observation per individual in pilot runs.These smaller datasets were created to compare the answers provided by the h-lik method against the standard (glmer) GLMM implementation in R. glmer uses Adaptive Gaussian Hermite (AGH) quadrature for the integration of random effects.This is a numerically accurate, but computationally intensive method that is used as a benchmark of alternative estimation methods for GLMMs [4,16].The number of nodes determines the order of the AGH scheme and can be increased to improve accuracy at the expense of speed.In these analyses we examined AGH quadrature with 1, 5 and 9 nodes.Note that the LA is formally equivalent to an AGH with a single quadrature node, so that the comparison between AGH1 and h-lik reflect the speedup due to the use of TMB.

Modeling the risk of dyskalemia via Poisson GLMM
The data   for these analyses were death events (  = 1) recorded in the Health Facts database.This outcome was modelled as a Poisson variable, by exploiting the link between Poisson models and analysis of survival time [1].
Features of interest that were coded as fixed effects were the K + level, the Charlson comorbidity index, a well validated instrument of the comorbid conditions affecting an individual, the prevailing level of renal function (estimated Glomerular Rate, eGFR, computed by the validated CKD-Epi formula), patient's age, type of healthcare facility (e.g.hospital, nursing home, outpatient clinic), patient's age, race and gender.Natural splines were used to probe non-linear relationships between the outcome of interest, renal function (eGFR) and the potassium level as well as age.We run two separate analyses with different random effects specifications: one in which the grouping level was the individual and one in which the grouping level was the facility.These analyses allow us to explore the scaling of performance against the dimensionality of the problem.The latter is largely determined by the number of random effects, since the vector of fixed effect coefficients is low dimensional, i.e. p=18 in all analyses.

Execution times of h-lik vs AGH based methods
For the timing information in the smaller datasets we run h-lik methods three times to generate timing information.We collected such information for all three steps, i.e. the generation of the tape, optimization and calculation of the Hessian at convergence.Timing information was generated in a high-end consumer workstation equipped with the i7-5960x octacore Intel Processor equipped with 32 GB of DDR4 RAM clocked at 2133MHz.All experiments were performed in Microsoft R Open v3.5.3.The base frequency of the i7-5960x is at 3GHz, but the processor was overclocked to 3.8 GHz for the timing experiments.In general execution time for these steps was not very variable (less than 15%).Therefore, only averages for these runs are reported here.Figure 1 shows the relevant timing comparisons.

Figure 1: GLMM execution times
Irrespective of the size of the dataset, the h-lik implementation was much faster than the AGH glmer implementation.All analyses with glmer generated warnings about at least one non-positive definitive Hessian matrices suggesting that the optimization algorithm got stuck close to the true global optimum.Furthermore, AGH1 failed to converge at all, when the measurements were grouped by the individual.In contrast to previous reports [5], we could obtain a solution from glmer for the 10% dataset for very large number of random effects but only if the AGH5 or AGH9 method were used.An Analysis of Variance (ANOVA) showed that the size of the dataset (F statistic 68.07) and the estimation method (F statistic 18.1) were much stronger predictors of these timing differences than the grouping structure of the data (F statistic 3.68).
Examination of the timing of the three steps of the h-lik method showed that the greatest amount of time was spent on optimization, e.g. for the 10% dataset with random effects at the individual patient level, tape generation took 137sec, optimization 3327 sec and Hessian calculation/uncertainty quantification 1,844 sec.
For the 1% dataset with random effects at the health facility level, the corresponding times were: 8 sec (tape generation), 17 sec (optimization), 1.7sec (uncertainty quantification).

Estimates of the-lik against AGH based methods
Convergence failure of the AGH1-9 methods implied that their estimates may not reliable and cannot be used as the sole basis of comparison for the h-lik estimates.We thus considered an alternative: the "zero" order LA that is also available from within glmer [16].The "zeroorder" LA is less accurate on theoretical grounds, since estimates will contain a definite, yet variable and unknown bias.Nevertheless, this approximation usually converge without warnings, and it did so in our datasets.The results of h-lik, AGH0, AGH9 are shown in Figure 2 as estimate (point) and 95% confidence interval for the two groupings and the two datasets.In general, the estimates of h-lik similar to the AGH methods, suggesting that the nonconvergence reported by the latter would not affect inferences about the fixed effects.The difference among the different methods are well within the uncertainty of the estimate; in fact, estimates were identical for most covariates up to two significant digits.
Nevertheless, there were substantial differences among the methods when the estimates of the variance components and the group specific random effects were examined.In these analyses, there is only a single variance component, i.e. the standard deviation of the random effect.Table 2 shows the estimates (and for h-lik the standard errors) for h-lik, AGH0 and AGH9.
Increasing the size of the dataset results in a greater agreement among the three methods and reduction in the magnitude of uncertainty in the estimates.Estimates by the AGH5 method were identical to the AGH9 in three significant digits (not shown).The distribution of the random effects was also similar among h-lik and the AGH methods (Figure 3).Similar results were obtained when the distribution of the random effects for the analysis that grouped observations by participant (not shown).As a final check of the calculations we repeated the analysis of the 1% dataset from a Bayesian perspective.This is justified because the system of equations ( 1) and (2) corresponds to a Bayesian analysis with noninformative (constant) priors for the fixed coefficients and the variance components.However, the Markov Chain Monte Carlo (MCMC) methods used in Bayesian analyses do not depend on the LA for the evaluation of the relevant integrals via the LA.Hence contingent on the convergence of the MCMC methods, Bayesian analyses may be considered the benchmark against which other methods of integration can be judged.For these comparisons we used the R package rstanarm.We only analyzed the 1% dataset at the facility grouping level and simulated four independent Markov chains in parallel.These simulations took 20,517 seconds to generate 1000 samples from each chain.The estimate of the variance component was 0.435, an estimate that was in rough agreement with the estimates generated by all the methods in Table 2 , but was numerically closer to the hlik estimate of 0.422.We also compared the random effects estimates from h-lik and the AGH methods against the median estimate of these quantitates from the MCMC (Figure 3) output.These analyses show that the AGH0 generated the large relative average difference (gray) from MCMC compared to the other methods, while h-lik 's performance was similar to that of the AGH9 method.

Figure 3. Random Effect h-lik/AGH estimates against MCMC
In summary: the h-lik method can fit very large GLMMs with accuracy that is equivalent if not better than the state of the art AGH methods.For the problems we examined, the h-lik gave results that are indistinguishable from the gold standard of MCMC integration but at a fraction of the execution time.

What is a healthy potassium level?
Having established that the h-lik implementation generates reasonable results in the smaller Health Facts datasets, we now turn to the analyses of the entire data.
For these analyses we used an Xeon E5-1603 quad core processor clocked at 2.8 GHz and equipped with 128 GB DDR4 RAM (clocked at 2133 MHz).The operating system for these analyses was Ubuntu 16.04 Memory was the major limiting resource for these analyses.Generating the computation tape consumed 238 GB (RAM + swap file).On the other hand, during optimization only 60GB of RAM was used.As a result of the heavy swap file use, timing, running the analyses on the entire dataset took nearly 36 hours.The major output was a response curve relating the relative risk of death based on a measured potassium level (Figure 4).This relationship was adjusted for a large number of individual characteristics to reduce confounding by other patient factors e.g. higher comorbidity or advanced age.In addition to these global features, we used random effects to capture residual heterogeneity at the patient (or the healthcare facility that evaluated the patient) level.To generate the graph, we used a reference potassium of 4.5 meq/l (this is the middle point of the reference range used by most clinical laboratories).The figure shows that risk of death numerically close to one over a broad range of potassium levels.Only outside the range of 3.5 -5.5 meq/l did the risk exceed 10%.This relationship was qualitatively (Ushaped) and quantitatively similar to the one previously reported in a much smaller study, reaffirming this finding and extending it throughout the broad range of patients seen in the diverse healthcare facilities using Cerner.Individual variability around this curve was small and quantitatively similar to the estimate based on the 10% dataset in Table 2.

Significance and Impact
In this short report we demonstrate feasibility of fitting GLMMs to mine large clinical EHR databases.This is an area in which many existing implementations of GLMM are found to be lacking.There are several implications for future work are noted below.

Random effects models and very big data
GLMMs are one of the most versatile modeling techniques in statistics.In recent years, it has been noted that numerous models used for the non-parametric, flexible modeling of data can be viewed as special cases of the GLMM.Such models include penalized spline approaches in statistics, but also "workhorses" of Machine Learning such as kernel based methods [19].However, the applications of GLMMs to big data to date has been limited by the computational bottlenecks inherent in evaluating the high dimensional integrals that arise in their mathematical formulation.A previous biomedically oriented paper reported difficulties in fitting random effects models with 890,934 random effects using the implementation of glmer in 2012.The authors resorted in a divide and conquer technique to fit meta-regression models to the same data.However, because of their approximative nature, any divide and conquer technique will be accompanied by a nonnegligible amount of bias.Furthermore, certain approximation approaches will exhibit a dramatic loss in statistical efficiency, discarding information and essentially converting a "big", information-rich dataset to a much smaller one [16].Having the ability to fit a large dataset with the same rigorous techniques that work in small dataset would thus represent a major advantage, because of the preservation of statistical efficiency.In our analyses we effortlessly fit a problem that was 16 times as large as the one previously reported in the biomedical literature.We feel that the slice of the Health Facts that generated for this project could be a useful benchmark problem given its challenging and unbalanced repeated measures structure for similar biomedical analyses.Although GLMMs are frequently encountered in the analyses of biomedical data, their scope is certainly not limited to this application domain.For example, an implementation of GLMMs by the LinkedIn software engineers showed that the random effect models performed among the top 5 to 10 methods in numerous nonmedical datasets.Taken together the results of the LinkedIn team and our analyses, suggest that GLMMs deserve a closer examination in this era of massive data.

Implementation of GLMMs
We adopted a computationally efficient implementation of GLMMs by mapping of the hierarchical likelihood approach to the software capabilities of the TMB library in R. The latter is based on a Maximum Likelihood (rather than REML) estimation procedure for GLMMs so that our mapping of theory to software was indeed necessary.TMB uses state of the art algorithmic differentiation methods and the Laplace approximation to implement Maximum Likelihood for these models.The package glmmTMB was recently reported [3] to be able to fit Random effects models for GLMMs using an interface similar to R's glmer.However, to our knowledge this package does not implement the h-lik approach, although it does appear to be able to estimate models via REML.We feel that a major advantage of the original interface of TMB is that one can flexibly swap optimizers depending on the size of the problem at hand.For example, LinkedIn's parallel optimization approach for GLMMs could be integrated with TMB's approach to mixed model estimation to develop extremely scalable implementations that execute across (possibly) heterogeneous computing infrastructures.Future work should also examine the implementation of GLMMs within more mainstream Machine Learning platforms such as TensorFlow.In fact, AD of machine learning models and scalable optimization is what TensorFlow excels at.Such implementations may have a more favorable memory footprint than the one in the TMB implementation, reducing the hardware requirements needed to analyze EHR data.Therefore, research in alternative implementations of the h-lik has the potential to make the theoretical soundness of GLMMs available to practicing data scientists in addition to statisticians.

But what did we learn about the potassium level?
This work was motivated by a clinical question that the kidney specialists' authors posed internally a few years ago.We feel that the analysis presented herein answered this question by re-affirming the "U" shaped curve between the potassium level and mortality previously suggested [13].On the basis of our data, potassium levels between 3.5 to 5.5 meq/l are unlikely to be associated with excess mortality.Interindividual variability around this relationship appears to small Hence, our bracketing of a globally valid reference range for potassium levels could be used to inform the design of clinical trials to guide the proper use of potassium supplements and potassium lowering drugs.

Figure 2 :
Figure 2: Estimates and 95% confidence intervals for the fixed effect components

Figure 4 .
Figure 4. Mortality across the range of potassium levels