Survival Analysis with Electronic Health Record Data: Experiments with Chronic Kidney Disease

Generalizability.

An EHR population is often tied to a specific location, via the so-called “catchment area” of the governing institution. The size and socio-demographic makeup of the catchment area depends on the institutional infrastructure, and economic as well as administrative issues. These issues may manifest themselves in practice as differences in distributions of risk factors, such as race, socio-economic status, and institutional and medical protocols, which could result in EHR-specific disease case severity and risk profile mixes. When not accounted for, these differences could lead to analyses and results for the specific EHR population being non-representative of the population at large [10, 11]. In addition, uncertainty in patient socio-economic predictors corresponds to the measurement error problem, which needs to be dealt with statistically as described in [11].

Simple examples of this include differences in EHR populations due to institution’s physical location, such as rural versus urban. Location may be a proxy for access to care, as well as environmental and socio-economic factors. More complex examples include situations where prevalence or other aspects of the disease under study are linked to the EHR population characteristics in a confounded way, where confounders are unmeasured and potentially time-varying [12, 13], such as relating health care patterns and the number of visits a subject has made to the doctor to disease severity. Many of these issues are not new, and have been extensively dealt with in the area of environmental health, health economics, and healthcare policy, especially outcomes research [14] and provider profiling literature [15, 16].

Issues such as these underline the importance of sophisticated statistical modeling as well as replicability for EHR at multiple institutions, and in multiple socio-economic strata. Development of formal methods for comparisons of specific EHR results to other EHRs as well as results published in the literature from alternative types of studies, perhaps using meta-analysis tools and Bayesian techniques [17, 18], could shed additional light on these complex issues.

Cohort definition.

EHR contains information on a variety of patients, differing in socio-economic and health histories, health seeking behavior, as well as types and quantity of data in their records. The act of performing an EHR-based study requires an algorithm-driven extraction of a specific and well-defined EHR subpopulation; including all patients from EHR is neither practical nor advisable. This process requires decisions regarding inclusion and exclusion criteria, length of record, and amount and quality of data in the record. The final cohort of patients is ideally well defined, and representative of the larger specific subpopulation it was chosen from.

Identifying patient cohorts from EHR data is an open research problem [19]. For some conditions, selection of patients can rely on simple laboratory test queries. For instance, patients with chronic kidney disease can be selected by their eGFR (estimated glomerular filtration rate; this variable is computed based on creatinine levels, a routine laboratory test) because eGFR is often the sole variable used to diagnose CKD. For other conditions, however, selecting patients is more complex, as the clues to the presence of disease are spread over heterogeneous data elements in the record (e.g., laboratory test values, mentions of symptoms in the notes, and presence of specific medications). Furthermore, an institution’s record may miss relevant information about a patient history. For instance, a particular medical event which should preclude a patient from being identified as a case might have not been recorded into the EHR because the patient sought care in a different institution. As such, there is often uncertainty associated with the definition of cases and control among the patients. A careful cohort definition aims to reduce this uncertainty, ultimately increasing the usability of the findings [11]. However, when not possible, misclassification in cohort definition can be statistically incorporated into effect estimates as described in [11].

Missing data.

EHR data are generally collected over long time periods, which means that missing data is one of the most important statistical issues to consider. Missingness in each of the three processes (disease related, patient related and healthcare related), has a potential for biasing the results of statistical analyses [11, 20, 21]. Additionally, uncertainty about the missingness mechanisms themselves presents another layer of complexity for statistical analysis of EHR data.

Amongst many types of missingness, record gaps (incomplete data) and censoring (left, right, and interval) are some of the most commonly encountered in EHR. More formally, if we let an individual’s true health status be represented by a multivariate continuous time stochastic process X(t), his or her interaction with a healthcare system as another stochastic process Y (t, X(t)), then the observed EHR record is a third stochastic process h(X(t), Y (t, X(t))). In this framework the first processes is latent (hidden, unobserved), while the second and third one are only partially and often irregularly observed.

In many scientific contexts, missingness might be relatively benign and non-informative, and manifested simply as non-uniformly sampled data in a way that is not related to the disease of interest (missing at random). However, in the EHR context, the missingness can be informative and, more specifically, related to the patient’s disease status itself. In those situations, postulating the mechanisms underlying the missingness, as with Bayesian approaches [22], and formally assessing robustness of results with respect to the implicitly defined missingness mechanisms, via subsampling and comparisons with other large EHR populations and clinical studies, is paramount.

Censoring.

Real patients move, change insurance, and are seen by different healthcare centers and institutions that collect data independently. Moreover, real patients vary in healthcare seeking behaviors, both reactively and proactively. This means that measurement pattern (visits as well as lack of visits) is sometimes related to health state. However, it is unknown whether lack of measurement implies absence of the disease, a severe case of an unrelated morbidity or mortality, a severe case of the disease itself treated by a different healthcare provider, or is a reflection of a patient’s attitude towards healthcare. In the time-to-event context, this means that censoring is a complex issue, where different types of censoring (left, right, interval censoring) may occur for each patient, and cannot always be assumed to be a separate, random process. For example, in the study of progression of chronic kidney disease, assuming onset of disease along with subsequent information are indeed present in a patient’s record, lack of a progression event causes the patient’s time to progression to be right censored at the last recorded visit date. Put differently, a lack of measurement does not imply health, and the information in the record only informs the lower bound on the time to progression.

Like in patient cohort selection, one way to alleviate the effect of this challenge is to focus on a patient population for which EHR has enough longitudinal and representative information. This may be an attainable goal for EHRs from institutions with captive populations which do not tend to move (long-term residents), or where patients are observed for all care, as may be the case with inner city populations and community hospitals.

Variable selection.

EHR data exist in the form of laboratory and other physiologic measurements, structured and unstructured notes, orders of various types (e.g., medication orders), billing codes, and images (e.g., CAT scans). Within each of these variable types lie hundreds to thousands of concepts and variable types. Beyond the sheer number of variables to choose from, not all variables are measured with the same frequency, on the same severity scale, or with the same quality and certainty. For example, an event is sometimes reported in a clinician’s note with qualifiers such as “perhaps mild anxiety” or “patient denies smoking”. These issues make the automatic variable selection process in the EHR setting remarkably complex.

Additionally, variables in the record may have different granularities and definitions. For instance, laboratory measurements for the same type of test, such as glucose, are often recorded under different measurement names depending on the healthcare setting where the measurement occurs. In variables extracted from notes, granularity issues are even more pervasive – a patient’s “cough” and “productive cough” could be considered the same concept in a study of kidney function, but should probably be considered separate variables when studying pneumonia [23]. Even further, some measurements become obsolete over different periods of time. For example, a troponin measurement is representative of a physiological state for a day, while a hbA1c value is valid for six months. The dimensions that capture the human condition, and the healthcare system itself, are continually evolving.

While there exist advanced methods for variable selection catered to large data (for example, see [24] for methods for variable selection for survival analysis), none developed to date are especially equipped to deal with EHR-specific issues. The variable selection step is a critical step to EHR-driven analysis, and research into adaptive and time-varying variable selection methods would be of great interest.

Statistical models: joint survival and longitudinal modeling.

A significant amount of EHR appeal lies in its ability to provide joint modeling of patient’s time-varying processes and disease progression. As the medical field moves towards personalized medicine, ability to provide patient-specific counseling based on all available patient’s data would be a powerful asset to any EHR based system. Joint modeling of longitudinal and time-to-event data has been an active area on the forefront of biostatistical research, with many notable contributions over the past two decades (see for example [25–28]).

Within a large EHR, considerable diversity can often be observed in the frequency, quality, and type of longitudinal measurements, as alluded to in the earlier missing data discussion. Many individual patients have a burst of measurements over a short period of time, while allowing years to pass without being seen by the provider. A fraction of patients have routine regular measurements, and aperiodic bursts of dense measurements; see Fig. 1 in [29] as an example. Yet another fraction has very sparse, irregular measurements, with no particular discernible pattern.

One way to deal with irregular measurements would be to follow the aggregation method in [30]. Alternatively, interpolation or latent variable modeling of covariate processes may provide a solution for imputing the missing measurements over time. There is a history of modified interpolation working well in the EHR context [31]. Nevertheless, interpolating many months or years of laboratory test values can lead to serious problems, not the least of which involve drastic changes in the derivatives of a patient’s health state. As an example, consider a patient whose blood glucose level is normal for two years, but is followed by a sudden increase at the end of the two year period. If the patient’s glucose is only measured at the beginning and the end of the two year period, linearly interpolating between the two measurements will obfuscate the sharp change likely masking (left-censoring) the true onset of a disease. Researchers have have tried to understand how to deal with this problem in the EHR context, using methods ranging from “half-life” analysis of laboratory measurements [32], to information entropy for establishing measurement rates such that they remain reliable [33]. Research on joint modeling for observed covariate processes and disease progression in the EHR context promises to be a fruitful area of research.

3.1. Chronic Kidney Disease

Chronic Kidney Disease (CKD) is a progressive condition in which kidney function decreases over time. When left untreated or poorly managed, CKD can lead to renal failure requiring dialysis and transplantation. In the United States its prevalence is reaching epidemic proportions (13% of the general population have a diagnosis of CKD, and 6.3% have stage 3 or higher CKD), and it costs the nation billions annually, totaling 17% of total Medicare expenditures for patients of age 65 and over [34].

The National Kidney Foundation has established five stages of progression of CKD to reflect levels of kidney damage and ensure that each stage gets appropriate care. Staging is based on the glomerular filtration rate (GFR), which is considered the best estimate of kidney function. In practice, the estimated GFR (eGFR), calculated using routinely collected variables, such as plasma creatinine level, age, and gender, is an acceptable metric for kidney function.

The official definition for CKD stage 3 (or moderate CKD) is kidney damage as reflected by a moderate decrease in eGFR to the range of 30 − 59ml/min/1.73m² for at least three months [35]. Patients are recommended to be referred to a nephrologist when reaching stage 3, so that an appropriate treatment can be devised to slow down the decrease in kidney function. Stage 4, or severe CKD, is defined as eGFR in the range of 15 − 29ml/min/1.73m². Stage 5 (GFR below 15ml/min/1.73m²) is considered kidney failure. If CKD is left untreated, eGFR decreases in a somewhat linear fashion [35]. In our study of CKD progression, we thus focus on the survival analysis of progression from stage 3 to stage 4.

Many conditions can lead to CKD, including diabetes, cardiovascular disease, and hypertension [35]. The mechanisms of damage vary – a number of factors have been found to be predictive of a decline in function and some factors may be protective [36–38]. Preventing dire outcomes requires both early detection, aggressive management plans, and resources [39]. Overall, CKD is a complex disease because decrease in kidney function affects many systems in the body, which in turn damage the kidney further. As such, progression of CKD is still poorly understood [38].

3.2. Dataset and Data Preprocessing

After IRB board approval, the dataset was assembled from a single institution, NewYork-Presbyterian (NYP) hospital, in New York City. Part of NYP is an outpatient clinic, which provides primary care. Patients who visit the clinic tend to be local to the hospital. Thus, in addition to routine visits to their primary care providers and consultation with specialists, their records are also more likely to contain information from visits to the emergency department and hospital admissions. Choosing a stable subpopulation like this one reduces censoring and informative missingness, and increases the likelihood of EHR containing most of the patient’s medical record. Furthermore, for the sake of a survival analysis study on a disease such as CKD, for which progression occurs at the resolution of years, our dataset is comprised of regular patients from our outpatient clinic.

While the EHR contains much data, we focused on three types of variables: demographic variables, healthcare variables, and laboratory test variables. Healthcare variables, such as number of visits, are particularly attractive in investigating EHR data, as their dynamics reflect a mixture of processes: the disease-related process (more severely sick patients are likely to visit their doctors more than less severe patients), the healthcare process (routine visits and consultation visits often follow schedules as suggested by clinical guidelines, but also depend on access to care and types of insurance benefits), and the patient process (some patients are more adherent to care than others, and thus are more likely to visit their doctors). The EHR is the perfect conduit for capturing visit patterns. We hypothesize that healthcare variables, as derived from the EHR, provide valuable insight into disease progression. As for laboratory variables, they have been the primary types of variables investigated in disease progression thus far and as such made much sense to be included [38]. We leave for future work the use of other EHR variables, such as diagnosis codes and medication orders.

4.1. Multi-Resolution Hazard Estimator

Bayesian multiresolution models have recently been adapted to hazard functions, and equipped to deal with censoring and truncation [67–69]. Following [66], these models assume a binary partition tree structure for time-axis partition. A Gamma prior is placed on the total cumulative hazard function over a finite time interval, and a set of Beta priors are used to describe the recursively defined partition parameters (also referred to as the “split” parameters). While essentially providing a piece-wise continuous prior on the hazard rate, the tree structure allows several interesting properties that are not common to other piece-wise hazard models. This includes the ability to control the smoothness of the resulting hazard rate via hyperparameters of the Gamma and Beta priors. The MRH model was used to model a breast cancer recurrence risk across multiple clinical centers in [68], prostate cancer [71], and multivariate hazard for different subpopulations in [69].

In MRH, we define the “final time resolution” comprising J equal length time intervals ${\left\{\left(t_{j-1},t_j\right)\right\}}_{j=1}^{j=J}$, covering the time span of the study (t₀, t_J). The final time, t_J, is set at, or slightly before, the last observed failure in the sample. The piece-wise approximation to the hazard rate is done via a set of hazard increments, ${\left\{d_j\right\}}_{j=1}^{j=J}$, where each d_j represents the aggregated hazard rate over the interval (t_j−1, t_j), i.e. $d_j={\int }_{t_{j-1}}^{t_j}h(s)ds\equiv H\left(t_j\right)-H\left(t_{j-1}\right)$. To facilitate the recursive binary partition, we let J = 2^M, where M is a positive integer describing the “height” of the partition tree. M can be selected via model selection criteria as in [68], or pre-set using clinical input, as in [70, 71]. From a statistical point of view, M is usually chosen so that each time interval has multiple event occurrences within it.

For simplicity, let H denote the total cumulative hazard H(t_J) over the study interval (t₀, t_J), and let H_M,0 ≡ d₁, $H_{M,1}\equiv d_2,\dots ,H_{M,2^M-1}\equiv d_J$. For each m, m = 1, 2, …, M −1, we recursively define the “level-m” hazard increments as H_m,p ≡ H_m+1,2p+H_m+1,2p+1, where p = 0, 1, 2, …, 2^m−1. Then the corresponding split variable R_m,p are defined as H_m,2p/H_m−1,p. From here, the MRH tree can be specified by H, and the set of splits $R_{1,0},\dots ,R_{M,2^{M-1}-1}$, and any d_j can be represented as a product of H and a certain branch of split variables. For example, when M = 3, we have d₁ = HR_1,0R_2,0R_3,0, d₂ = HR_1,0R_2,0(1−R_3,0), …, d₈ = H(1 − R_1,0)(1 − R_2,1)(1 − R_3,3).

We adopt Beta priors for split parameters R_m,p’s, and a Gamma prior on H as in [68]. For example, the priors for H and R_m,p in a M = 3 (J = 8) model would have the following priors:

With this parametrization, the prior expectation of each hazard increment is determined by the shape parameters of the Beta priors, and the hyperparameters of H, as each increment d_j is a product of independent variables; in the M = 3 example, E(d₁) = E(H)E(R_1,0)E(R_2,0)E(R_3,0). The tree hyperparameters are thus directly linked to the distribution of the hazard increment, specifying a-priori expected values of each hazard increment as well as the correlation between the hazard increments [67, 71], thus controlling the smoothness in the multiresolution prior.

4.2. Estimation

The multiresolution prior for the hazard rate can be used as a prior for the baseline hazard function in any model for time-to-event data. One of the most commonly used models in survival analysis is the Cox proportional hazard model [45], where the hazard rate for patient i is modeled as follows:

Here, the hazard rate for a specific observation i is obtained via multiplication of the baseline hazard rate h₀ (the hazard in the baseline group of patients) by $\text{exp}\left(X_i^{\prime }\beta \right)$. The parameter vector β is known as log hazard ratio. The log hazard ratios are the effects of covariates in X_i, a vector of p explanatory variables X_i,1, X_i,2 ···X_i,p for observation i. The hazards in different patient groups are thus assumed to be proportional to each other over the entire study period, giving this model its name.

Unlike the Cox proportional hazards model [45] which aims to estimate the effects of covariates while treating the hazard function as a nuisance parameter, the MRH model estimates the hazard and all parameters jointly, allowing for joint inference which is based on the joint posterior distribution estimated using MCMC. The likelihood function of the MRH proportional hazards model based on right censored data from N patients is:

Here, δ_i = 0 if T_i is right censored, and δ_i = 1 if not (i.e. if that patient did not have an observed event within the study period).

The Bayesian MRH model is then completed by specifying prior distributions for β, as well as on hyperpriors in the MRH tree parameters (a, λ, k, and the set of γ_m,p parameters). In what follows, a flat prior shape on the baseline hazard function is assumed, via setting k = 0.5, and γ_m,p = 0.5 for all m, p, making the prior on the split parameters a symmetric beta density, R_m,p Beta(a/(2^m), a/(2^m)), evenly splitting the hazard among the intervals. The Gibbs sampler [72] steps for the remaining variables (H, {R_m,p}, a, λ, and β) are the same as those in [68]. We briefly outline each full conditional distribution below, making up the core of the Gibbs sampler. In what follows, the superscript ‘⁻’ is used to denote a vector of all other parameters and data except for that parameter:

1. If a is given a vague zero-truncated Poisson prior, $\frac{e^{-{\mu }_a}{\mu }_a^a}{a!\left(1-e^{-{\mu }_a}\right)}$, the full conditional distribution for a is given as follows:

   $$\pi \left(a\mid a^{-}\right)\propto \frac{1}{{\lambda }^a\Gamma (a)}{H}^a\frac{{\mu }_a^a}{a!}{\Pi }_{m=1}^M{\Pi }_{p=0}^{2^{m-1}-1}\left\{\frac{R_{m,p}^{a/2^m}{\left(1-R_{m,p}\right)}^{a/2^m}}{\text{B}\left(a/2^{m}a,a/2^m\right)}\right\},$$

   where B(a/2^ma, a/2^m) is the Beta function.
2. If the scale parameter λ from Gamma prior for the cumulative hazard function H is given an exponential prior with mean μ_λ, the resulting full conditional is:

   $$\pi \left(\lambda \mid {\lambda }^{-}\right)\propto \frac{1}{{\lambda }^a}{e}^{-\frac{H}{\lambda }}{e}^{-\frac{\lambda }{{\mu }_{\lambda }}}$$

3. A normal $N\left(0,{\sigma }_{\beta }^2\right)$ prior on each log hazard ratio, β_j, leads to the the following full conditional distribution:

   $$\pi \left({\beta }_j\mid {\beta }_j^{-}\right)\propto \text{exp}\left\{-\frac{{\beta }_j^2}{2{\sigma }_{\beta }^2}\right\}{\Pi }_{i=1}^N\left\{{\left[\text{exp}\left(X_{i,j}{\beta }_j\right)\right]}^{{\delta }_i}\text{exp}\left(-\text{exp}\left(X_i^{\prime }\beta \right)HF\left(T_i\right)\right)\right\},$$

   where F(T_i) = H(min(T_i, t_J))/H(t_J).

4. The full conditional for H, π(H|H⁻) is a Gamma density, with the mean of $\frac{a+{\sum }_{i=1}^N{\delta }_i}{\left[(1/\lambda )+{\sum }_{i=1}^N\text{exp}\left(X_i^{\prime }\beta \right)F\left(T_i\right)\right]}$, and variance of $\frac{a+{\sum }_{i=1}^N{\delta }_i}{{\left[(1/\lambda )+{\sum }_{i=1}^N\text{exp}\left(X_i^{\prime }\beta \right)F\left(T_i\right)\right]}^2}$.

5. The full conditionals for each R_m,p are given as follows:

   $$\pi \left(R_{m,p}\mid R_{m,p}^{-}\right)\propto {\Pi }_{i=1}^N\left\{{\left[h_0\left(T_i\right)\right]}^{{\delta }_i}\text{exp}\left(-\text{exp}\left(X_i^{\prime }\beta \right)H\mathcal{F}\left(T_i\right)\right)\right\}{R}_{m,p}^{2\gamma k^{m}a-1}{\left(1-R_{m,p}\right)}^{2(1-\gamma )k^{m}a-1}.$$

The sampling from full conditionals is done either directly, or via a Metropolis-Hastings like algorithms. All estimation was performed using a custom made MRH library in R [73]. For each model we consider, statistical inference was based on 5 chains, each with 200,000 iterations, the first 50,000 iterations burned, and the remaining iterations thinned to remove autocorrelation. Convergence was determined via trace plots, autocorrelation plots, and the Geweke test statistic [74]. Parameter estimates were calculated using the median of the posterior distribution of each parameter, generated by the Gibbs sampler chains. Uncertainty was communicated via 95% central credible intervals, whose bounds were calculated as the 2.5 and 97.5 percentiles of the simulated posterior distribution.

5.1. Experimental Setup

The goal of the analysis is to identify variables associated with progression of CKD, as derived from our patient cohort. Because of the nature of the EHR dataset, we experimented with different settings for the survival analysis and survival models. Table 3 shows descriptive statistics of the different settings.

5.2. Basic Healthcare Model

Differences between the two definitions of the observation time were minimal, based on the summaries in Table 3 and the MRH model results in Table 4. Summaries of the difference between the lower and upper bound definitions show an average difference of less than 5 months for all groups, with a maximum difference of 3.5 years. In addition, the MRH model results show very minor differences in the parameter estimates of all models, particularly in the comparison of the 95% credible intervals. It does not appear that the different definitions of the observation time have a large impact on the model estimates. For this reason, moving forward our discussion of the results will only focus on the models using the upper bound definition of the observation time.

Figure 1 shows the smoothed estimated hazard rate for the three different patient subgroups, using the second definition of the number of visits (models A2, B2, and C2). As expected, the hazard rate for all patients is higher throughout the course of the study, and especially so within the first few years since stage 3 diagnosis (as recorded in the EHR), when compared to the 3m and 6m groups (patients who had at least one laboratory test measure at least 3 (6) months prior to stage 3 CKD diagnosis). However, there appears to be almost no difference in the hazard rate between the 3m and 6m subpopulations of patients. For all groups, we observe that the hazard rate is highest in the first year, then flat between years 5 and 10 post stage 3 diagnosis, followed by a slight increase after ten years.

Figure 2 provides the summary of the results in terms of probability of developing stage 4 CKD within one year from stage 3 CKD diagnosis. Among all patients this probability is centered around 3.0%, which is approximately twice of that for patients with at least one laboratory test measure recorded in the EHR at least 3 (6) months before the diagnosis date. The probability in these two groups, 3m and 6m, is centered around 1.5%, confirming once again that there is little practical difference among 3m and 6m patients. This figure also shows the uncertainty about these probability estimates, indicating that there is almost complete overlap in information about the probability of failure within first year from stage 3 diagnosis for 3m and 6pm patients, while all patients are clearly at an increased risk in spite of some distributional overlap with the other two groups.

With respect to the effects of age on the hazard rate, among all patients the older ones are found to have higher risk of stage 4 CKD. However, the opposite seems true when 3m or 6m patients are analyzed. For example, Figure 3 shows that in the model with all patients, probability of a 40 year old developing stage 4 CKD within 5 years since stage 3 diagnosis is centered around 5.3%, while this probability is centered around 6.3% for an 80 year old. However, for 3m patients, the probability of a 40 year old developing stage 4 CKD within 5 years is centered around 4.8%, and around 2.6% for an 80 year old. In general, censoring rates will have a dramatic effect on covariate effects and survival times [21, 75]. In this particular CKD case, there are two censoring issues: possible informative left-censoring for patients who were in EHR for less than 3 months, and right censoring of either older patients or those patients who were in the EHR a long time with censoring due to other diseases or death. Thus, one possible explanation for the estimated probability of stage 4 for older patients in the “All” group might be the more rapid stage 4 onset, implying that the effects of age are not as impacted by censoring as they are in the 3m patient group.

In terms of the effect of number of visits, there were some notable differences between all patients and those patients who were in EHR for at least 3 months. For all patients, each additional visit to the doctor in the year prior to stage 3 diagnosis decreased the risk of stage 4 CKD by 11.3% (definition 1 for number of visits) or 13.9% (definition 2 for number of visits). However, the impact of the number of visits was more moderate for the patients who were in the EHR system for more than 3 months: for them, each additional visit to the doctor in the year prior to stage 3 CKD was associated with an increased risk of stage 4 CKD by 2% (definition 1), or decreased the risk by 7% (definition 2). This pattern can be observed in figure 4, which shows that the probability of developing stage 4 CDK within 5 years of stage 3 CKD diagnosis is highest for all subjects who had zero visits to the doctor in the year prior to diagnosis (using definition 1). However, this probability is lowest for the subjects who were in the system at least 3 months prior to diagnosis, and who did not visit the doctor (using definition 1). It is important to note that all the parameter effects’ 95% credible intervals included the null effect. All non-exponentiated covariate effect (log-hazard ratio) estimates, and the corresponding credible intervals, are shown in Table 4.

In terms of the effect of the number of lab tests, an increase in the number of tests (which can be a proxy for how sick an individual is) was associated with an increase in the hazard of stage 4 CKD. This, combined with the effects discussed in the above paragraph, indicate that while a large number of laboratory tests are associated with increased hazard of stage 4, frequent visits to the doctor can decrease the hazard of 4 CKD.

Finally, in terms of gender effects, male gender was estimated to be positively associated with the hazard of stage 4 diagnosis, across all models, although the only 95% credible intervals which did not include the null effect were obtained in models with all patients.

5.3. Clinical Model

Results from the basic healthcare models allowed us to conclude that the 3m subgroup is the most appropriate set of patients to use in studying the hazard of stage 4 CKD, as it seems that the data set contains a separate subgroup of patients who do not have a recorded true stage 3 diagnosis date. Additionally, as the laboratory tests are often performed over the course of a several successive days, it seems also reasonable to adopt the second definition of the number of visits. Starting from there, we extended this basic model (B2) to include 17 laboratory test variables defined in section 3, as well as the 17 indicator variables defined in section 3. The laboratory test variables considered (25-OH vitamin D, bicarbonate, blood urea nitrogen, calcium, ionized calcium, chloride, high-sensitivity C-reactive protein, hematocrit, hemoglobin, potassium, magnesium, microalbumin creatinine ratio, sodium phosphorus, parathyroid hormone, triglycerides, and uric acid).

To determine which of the 34 total laboratory variables (17 laboratory test variables and 17 laboratory test indicators) were most important to stage 4 progression, we performed a backward-forward Cox proportional hazards regression, adjusted for mandatory inclusion of non-laboratory variables already in the B2 model. The subset of variables selected this way (blood urea nitrogen, calcium, hematocrit, magnesium, microalbumin creatinine ratio, and sodium as continuous laboratory test covariates, and chloride, and C-reactive protein as indicators), along with gender, age, number of visits (type 2) in the year before stage 3 diagnosis, and the number of laboratory tests in the year prior to stage 3 diagnosis, were all included as predictors in the MRH model, which we refer to as “Optimal model”. We also consider a second model, which includes all laboratory indicators in addition to the selected laboratory variables and basic healthcare variables, as main effects. This second model is referred to as the “Full model”.

Figure 5 shows the hazard rate results from the MRH model: the hazard rate for the baseline group (females, 62.2 years old, with no visits or laboratory tests in the year prior to stage 3 CKD diagnosis) was estimated to be much flatter and lower in magnitude than the one shown for the B2 model. It is interesting to contrast these two models, as the optimal model is in fact model B2 with laboratory tests added as covariates. The baseline hazard rates for the two models look fairly different (Figure 5): the hazard rate of model B2 is lower than that of the optimal model, and the hazard rate of the optimal model shows an increase in the hazard rate after ten years which is not as visible in the hazard rate for model B2. Both models, however, suggest that baseline patients are at low risk for progression to stage 4. The hazard rate is again estimated to be higher for males, but the estimate of the effect is much smaller in this model when compared to the basic B2 model (Table 5). This is likely due to the fact that gender actually carries some of the same metabolic information as the laboratory test measurements, and that the laboratory test values themselves account for some of the differences between males and females. As with the basic healthcare model, the number of visits is estimated to have a negative effect on the risk of developing stage 4 CKD, when adjusted for laboratory variables. However, unlike in B2, in this model the number of laboratory tests is also estimated to decrease this risk, when adjusted for laboratory variables.

The optimal and full models show very similar estimates of effects for common variables. However, the effect interpretation in both models is slightly unconventional. The lab test value effects summarized under “Interaction terms” in Table 5 under the “Full model” column heading should be interpreted as the log hazard ratios of these lab test values conditional on the fact that the patient had that particular lab test ordered. For example, the full model shows that estimate of the indicator for hematocrit is equal to 0.251, and the estimate for the lab interaction term is equal to −0.358. Thus, for patients who had hematocrit test done in the year prior to stage 3 diagnosis, the estimated hazard rate is 29% higher than for patients who did not have any hematocrit measures ordered during that year, all else being equal. However, among those patients who did have a hematocrit test in the previous year, larger hematocrit values led to decreased risk (an estimated 30% decrease for every standard deviation increase on the log scale) when compared to patients who had lower hematocrit values.

Results from both lab models (Table 5) show that the largest impact on the hazard rate comes from increases in blood urea nitrogen and microalbumin creatinine ratio values. For example, patients who are 1 standard deviation above the mean for each lab (on the log scale) have a hazard rate more than 3 times higher than patients who had average measures for these two labs. This finding is not surprising as microalbumin creatinine ratio is well recognized as the first predictor of progression in the literature [76], and is one of the few labs which reflects kidney function directly rather than symptoms and quantities the kidneys help regulate.

Additionally, patients who had at least one high sensitivity C-reactive protein laboratory measurement during the year prior to stage 3 diagnosis have the hazard rate almost 4.5 times higher than patients who did not have that measurement during the year prior to stage 3. This finding too is not unexpected, as C reactive protein is a marker for inflammation and is typically ordered for patients with heart disease or chronic inflammatory conditions such as lupus, often seen as comorbidities of CKD. In essence, CRP could be an important proxy for a the presence of such comorbidities associated with elevated risk of CKD progression [77].

Furthermore, lower levels of hematocrit, magnesium, and calcium are found associated with increased hazard of stage 4 CKD in both lab models. In addition, both models also estimate that higher levels of sodium are associated with increased hazard of stage 4 CKD, which is not unexpected given that high sodium concentration is generally a proxy for poor hydration. However, there was little evidence that the effects of these two variables were different from 0, as indicated by the fact that both of their 95% credible intervals contained 0, for both the optimal and full model.

5.4. Comparison to Traditional Survival Analysis

So far we have discussed two types of MRH models: six basic healthcare models and two clinical models. While MRH is a powerful class of models that can be used to provide joint estimates of the hazard function and covariate effects, along with their uncertainty, a simple parametric model for the hazard function might still provide a reasonable alternative in some datasets. To test whether that is the case for the different types of healthcare models we consider in this paper, we utilize two parametric models models, assuming an Exponential and a Weibull distribution, respectively. The hazard rate is flat under the exponential distribution assumption for progression to stage 4 times, which resembles the estimated hazard shape starting with 18 months post stage 3 CKD diagnosis (Figure 5). The Weibull model is similar to the exponential model, but allows for more flexibility in the shape of the hazard rate. Parametric models were calculated using the survreg package in R [73]. The models were compared using Bayesian Information Criterion (BIC [78]), and a custom designed goodness-of-fit (GOF) measure as defined below:

where | · | stands for absolute value, I_i{advance to stage 4 CKD > t} is an indicator variable which equals 1 if the subject i advances to stage 4 CKD after time t post stage 3 diagnosis and equals 0 otherwise, and P(subject i advances to stage 4 CKD > t) is the model-based probability of the subject i advancing to stage 4 after time t post stage 3 diagnosis. This probability is calculated based on the estimated MRH parameters and covariate effects for subject i. Patients who were censored before time t were not included in the GOF calculation at time t. Therefore, n_t represents the total number of patients in the cohort minus the number of patients censored before time t. This way, the maximum value the GOF statistic can take is 1, and that maximum can only be achieved at the baseline time t = 0. Lower GOF values indicate a better fitting model.

Results from the model comparison (Table 6) show that the predictive abilities of the Weibull model and the MRH model are very similar across all time points, with the MRH models performing slightly better, although the Weibull model had lower BIC values. The lower BIC values for Weibull model were not unexpected, given that the MRH model estimated a fairly flat hazard function, and that the flexibility of MRH models requires them to have a larger number of parameters than Weibull. The exponential model performed worse than the other two models, with very high GOF values. While the MRH model performs the best in predicting the advance to stage 4 CKD, it seems that a Weibull parametric model is a reasonable alternative for estimation of the hazard rate from stage 3 to stage 4 CKD. The Weibull model is relatively easy to implement and communicate in an applied medical research. The exponential model does not perform very well in this analysis. In general, we do not recommend using a flat hazard in analysis of long-time studies, such as those based on EHR or longer clinical trials [70, 71], where the hazard function is likely to change over time due to changes in factors such as population structure, health care protocols, and prevalence of comorbidities.

The GOF table also allows us to compare the MRH models to each other. Not surprisingly, the model that performs best is the healthcare model with the laboratory tests. Despite the increased number of parameters, this model has the lowest BIC of all the MRH models, and performs the best in predicting time to stage 4 CKD. Because this model uses more patient-level information, it is better able to predict when each person will advance to stage 4 CKD. It is also not unexpected that the model with all patients is harder to predict accurately when compared to the models using 3m patients, as its GOF and BIC values are higher.

In addition to comparing the hazard rate estimates, we also observed that the covariate estimates and among the different MRH models, the Weibull model, and the Cox proportional hazards model were very similar, including the 95% uncertainty bounds.

Generalizability.

EHR data are connected to their physical environment, making generalizability an issue. The Columbia University Medical Center AIM clinic EHR treats patients in lower socio-economic strata living in northern Manhattan. There is a large amount of hispanic patients, for which CKD progression is poorly understood [79]. Our results could be specific to ethnicity, demographics and their cultural implications, and general environment such as air quality, housing quality, etc. One way to increase the generalizablity is through cohort selection. To further increase generalizability of our results, we carefully defined our cohort by narrowing the population by comorbidity. For instance, we removed transplant and HIV-AIDS positive subjects because these patients have kidney failure for known reasons and because these patients progression, induced by treatment of a particular disease, differs from canonical CKD patients. Similarly, we removed patients with acute kidney disease whose kidney failure occurs on a considerably different time scale than CKD. It is worth mentioning that CKD is a disease with a particularly direct method for diagnosis because eGFR, used to diagnose CKD, is easily calculated using routine laboratory tests. Not all diseases are this easily defined (e.g., congestive heart failure, pneumonia).

Data missingness.

When analyzing the longitudinal records of patients, missing data can occur because a patient seeks care in a different institution or because the patient does not seek care at regular times. To account for this, we selected a cohort which is local and comes somewhat regularly for primary care and acute care to our institution. Furthermore, the healthcare variables encapsulated the patterns of data missingness: number of visits made in the year prior to CKD onset, for instance, reflect the measurement patterns as well as potentially the patients’ attitude towards healthcare. While visitation rate does not fully account for the complex relationship between lifelong healthcare patterns and disease severity, the histograms of the visit data, stratified by length of the record, indicate that changes in patient visitation is a proxy for patient acuity. For example individuals with shorter record spans with a higher average numbers of yearly laboratory tests and visits (figure not shown) are more likely to become ill. Furthermore, the interaction laboratory terms allow us to determine how the hazard rate changes based on a patient having laboratory measurement at all.

Another type of missing data comes from the way information is recorded in the EHR. For instance, while race is a potentially important variable for CKD progression, we were not able to rely on it in the analysis, since race is not recorded for most of the patients in our cohort. This is a clear limitation to EHR data. In our current work, we are exploring ways to infer missing information from additional data source in the EHR, such as the clinical narrative about the patient, where race is mentioned.

Selection of variables that best predict outcomes.

Our variable selection procedure was governed by three processes: (i) biological rationale (to determine the 17 laboratory tests most relevant to CKD progression); (ii) logical thinking (to create a proper variable to denote the number of visits to the doctor); and (iii) basic statistical tools used to determine the optimal model. The results of this process allowed us to determine the effects of gender, age, and the number of visits to the doctor and number of laboratory tests in the year prior to CKD stage 3 onset on the time to stage 4. While forward-backward stepwise regression is a simplistic statistical tool for variable selection, we were able to identify laboratory tests that significantly effect time to stage 4 CKD, with results that match previously published clinical findings in the CKD medical community. In the future, we want to development and apply statistical tools for variable selection that would allow us to determine optimal models over a wider range of variables, including a larger set of laboratory tests, comorbidities, medications, and other potentially relevant patient features.

Censoring.

Because individuals are measured either at routine appointment times or when seeking care for an illness, it can be difficult to identify the exact start time of a disease — the time of clinical diagnosis can be days to months after the actual onset of the disease. Thus, left-censoring had an effect on the shape of the hazard rate. We accounted for these left-censoring effects by contrasting different subsets of patients selected based on the length of their observation record. The different subgroups did have different hazard rates, which implies that some of the subjects do not have their true stage 3 CKD diagnosis date recorded. In the future, we want to develop and apply methods for addressing uncertainty due to left-censoring in order to more accurately align patients or account for the varying uncertainty of diagnosis time for different patients, within a statistical model. We would also like to better infer, impute, or estimate the accurate time of diagnosis in particular.

Model selection and stability.

The MRH method is a robust technique for estimating the hazard rate from stage 3 to stage 4 CKD. Moreover the MRH was flexible enough to accommodate the complex fluctuations in the hazard rate while simultaneously estimating covariate effects. The data were best represented by the type III generalized extreme value distribution, the parametric Weibull model. In contrast, the exponential model was not a suitable representation of the data. This is an important point because it highlights the fact that the statistical model distribution used to represent variables collected in the EHR context must be chosen carefully and consciously. Even when a variable distribution could be governed by a known physiologic process, the health care process can intervene though selective measurement, netting a non-physiologic distributional representation of what should have been a physiologic variable.

Finally, the MRH results were consistent with the other modeling techniques, leading to the conclusion that the results do not depend strongly on the model. While the MRH was the best performer the Cox regressions revealed much of what was revealed though the MRH analysis. The MRH has potential to be much more representative, and much more able to account for EHR complexities. Nevertheless, when more simple modeling technology is desired, the Cox regression remains a good choice.

In the future, related specifically to MRH, there is much to be done to properly analyze the process of going from stage 3 to stage 4 CKD. Development of a model that can properly account for the uncertainty in left-censoring and interval-censoring of subjects is paramount to analyzing the hazard rate using all subjects available in EHR data sets. Use of all subjects may also differently inform us on how laboratory tests and other covariates impact time to stage 4 CKD.

A main focus of future work from the MRH modeling standpoint includes accounting for the longitudinal aspect of the laboratory measurements, both before and after stage 3 CKD diagnosis. This would allow clinicians to check the progress of the disease once a patient has been diagnosed with stage 3 CKD, and would allow them to sequentially update the estimates of risk of developing stage 4 CKD as each new laboratory test results arrives into the record. Such a methodological extension would allow for a more specific patient risk assessment, in line of the recent personalized medicine efforts. The MRH model can be extended to accommodate these types of measures, and will account for both the level of laboratory measurements as well as changes in the laboratory measurements, allowing researchers to identify higher and lower periods of risk based on gradual or sharp changes in clinical measures.