Large-scale collections of electronic health records (EHRs) are becoming useful for understanding disease progress, early diagnosis, and personalized treatments for many clinical diseases (Murdoch and Detsky, 2013; Hripcsak and Albers, 2013; Ghassemi et al., 2015a). EHRs contain rich patient information—disease history, demographics, vital signs, and lab results—that clinicians use to diagnose and treat patients. In this work, we are interested in developing a statistical framework that leverages medical data from a set of reference patients to enable personalized, real-time monitoring of new hospital patients. In particular, we consider data from the Hospitals at the University of Pennsylvania (HUP) containing hospital information for over patients, and the public Multiparameter Intelligent Monitoring in Intensive Care (MIMIC-III) data set with more than admissions from patients in intensive care units (ICUs) (Johnson et al., 2016).
One motivation for monitoring new patients is to characterize patient state to allow the early diagnosis of sepsis or septic shock. Sepsis is one of the leading causes of death in critically ill patients in the United States (Hotchkiss and Karl, 2003). Each year an estimated 750,000 cases of sepsis or septic shock occur in the US. The mortality rate of septic patients ranges from 20% to 30%, and accounts for roughly 9.3% of all US deaths (Angus et al., 2001; Kumar et al., 2011). Sepsis is usually developed during a patient’s stay in the hospital. However, accurate diagnosis of sepsis is difficult due to heterogeneous symptoms across patients (Pierrakos and Vincent, 2010).
One way to reduce the mortality rate of sepsis is to increase the accuracy of early diagnosis of sepsis. To do this, we might develop a model of patient state and fit this model to EHR data from previous hospital patients with sepsis. However, existing EHR data pose several challenges because they have been collected with traditional monitoring methods. Many of the covariates, lab results in particular, are sparsely sampled across patients. That is, there are only a small number of observations of any lab result per patient. For example, vital signs are generally taken once every three to four hours for inpatient data, and once every hour for patients in the intensive care unit (ICU). In contrast, blood tests requiring a blood draw are generally performed at most once a day. We see this sparsity in an example of 24 clinical covariates (Table 1) measured across time for a single patient, including four densely sampled vital signs (respiration rate, heart rate, systolic blood pressure, and body temperature) and 20 sparsely sampled lab covariates (Figure 1). Data missingness is systematic and not at random (Newgard and Lewis, 2015): a doctor will only order a test that will be informative in characterizing patient state relevant to diagnosis. Moreover, these time series data are not aligned across patients to a reference time point or disease onset; instead, patient intake is at time and release is hours or days later. The sparsity over patients and uncalibrated time series make the physiological progression of disease within patients or joint analysis of time series across patients challenging due to substantial uncertainty of patient state and rate of disease progression at any time.
|Vital||Respiration rate (RR)||87,076||493,964||147,445||291,466|
|Vital||Heart rate (HR)||96,317||527,989||227,951||294,746|
|Vital||Systolic blood pressure (SBP)||84,909||447,666||104,129||124,587|
|Vital||Body temperature (Temp)||80,597||364,286||94,468||56,533|
|Lab||Blood urea nitrogen (BUN)||12,528||71,825||21,751||25,102|
|Lab||Carbon dioxide (CO)||12,672||72,784||21,844||20,979|
|Lab||Glucose point-of-care (Glucose POC)||20,444||170,872||54,239||24,196|
|Lab||Mean cell hemoglobin (MCH)||12,587||69,736||18,379||20,877|
|Lab||Mean cell hemoglobin concentration (MCHC)||12,577||69,682||18,359||20,885|
|Lab||Mean cell volume (MCV)||12,587||69,751||18,380||20,875|
|Lab||International normalization ratio (INR)||5,733||38,810||17,005||15,735|
|Lab||Prothrombin time (PT)||5,722||38,844||17,007||15,734|
|Lab||Partial thromboplastin time (PTT)||5,872||41,894||19,596||17,185|
|Lab||Red blood cell (RBC)||12,600||69,776||18,387||20,876|
|Lab||Red cell distribution width (RDW)||12,580||69,757||18,381||20,877|
|Lab||White blood cell (WBC)||12,581||69,950||18,384||20,960|
In this work, we build a statistical framework that uses sparse, heterogeneous EHR time series data to monitor and predict vital signs and lab results for each patient in an online way. To do this, we first designed a nonparametric model based on Gaussian process (GP) multivariate regression to explore the correlations both within each clinical covariate across time and across clinical covariates given rich EHR reference data. Our model includes a highly structured GP kernel regularized using sparsity-inducing priors to avoid overfitting, allow interpretability, and ensure computational tractability. Second, we propose a framework based on nonparametric density estimation to tailor the empirical model to a patient-specific model for each new patient. For real-time monitoring, we update the empirical distribution from reference patients with patient-specific observations as measurements are observed. We evaluate our method, MedGP, on over 6,000 patients from three disease groups with more than four million measurements from the HUP data and one disease group from the MIMIC-III data set. We compare results to state-of-the-art approaches for patient online monitoring and investigate differences in correlations among covariates across disease groups.
2 Related Work
Related work falls into three areas of medical time series analysis: (i) incorporating noisy, heterogeneous, irregular, and sparsely sampled time series data; (ii) combining information across multiple time series; and (iii) exploiting reference data in addition to observations about the current patient to enable patient-specific predictions for a new hospital patient.
Most prior work has focused on modeling each clinical covariate separately. Due to the irregularity and temporal sparsity of medical data, conventional time series models, such as hidden Markov models (HMMs), autoregressive (AR) models, state-space models, and linear dynamical systems (LDS), are challenging to apply because of the assumption of regular measurement sampling in time. Recent work has focused on developing methods to compensate for the missing data in order to work with models that assume complete data. InKim et al. (2010)
-nearest neighbor (KNN) clustering were used for missing data imputation, and improvements in septic shock prediction were reported(Ho et al., 2014). In other work, a hierarchical switching LDS model was used to monitor the physiological signals during neonatal sepsis; the model allows the latent state of a patient to change during periods with fewer observations (Stanculescu et al., 2014)
. In an alternative approach, noisy and sparse time series data were smoothed temporally by putting Gaussian priors on the mean parameters of the Gaussian mixture model, which is related to a Gaussian process prior, although the distribution is over a finite-dimensional vector(Marlin et al., 2012).
Gaussian processes (GPs) are useful approaches for time series analysis because they can naturally capture irregular time series observations and estimate prediction uncertainties in a probabilistic framework (Roberts et al., 2012). For these reasons, GPs have been applied to the analysis of medical time series data. Previous work used a single-output GP regression model to smooth and impute each covariate independently (Stegle et al., 2008; Lasko et al., 2013). The Probabilistic Subtyping Model (PSM) added patient-specific information for smoothing temporal trajectories of clinical covariates and clustering disease subtypes (Schulam et al., 2015). PSM learns a mixture model based on a B-spline and GPs to impute the clinical measurements for patients with scleroderma. Demographic covariates, including gender, ethnicity, and clinical history, were also incorporated in the model. In an extension of PSM, the authors adapted patient-specific information to forecast specific clinical covariates (Schulam and Saria, 2015); the time series for each covariate was still modeled independently.
The idea of capturing the joint dynamics between vital signs and lab tests has also been explored. Using high-frequency regularly sampled time series, the dynamics between heart rate (HR) and blood pressure (BP) were modeled using a mixture of an LDS model (Nemati et al., 2012)
and a switching vector autoregressive model (SVAR)(Lehman et al., 2015). The joint dynamics estimated across covariates were reported to be associated with hospital mortality. In other work (Rizopoulos and Ghosh, 2011), a multivariate spline-based approach with linear mixed effects was used to predict multiple longitudinal outcomes and time-to-death of patients. Time series graphical models (TGMs) (Dahlhaus, 2000; Tank et al., 2015) have also been studied and applied for analyzing multivariate medical time series of ICU patients (Gather et al., 2002). TGMs model the partial correlations between each dimension of the multivariate time series as an undirected graph. However, both TGMs and SVAR models follow the assumptions of vector autoregressive (VAR) models, and thus assume the sampling interval of the time series is fixed across dimensions. In practice, this means missing data imputation needs to be done in advance (Tank et al., 2015).
Several multi-output GP frameworks have been proposed for other application areas. In the geostatistics literature, the linear model of coregionalization (LMC) characterizes correlations between outputs through a set of kernels and coregionalization matrices that estimate weights for pairwise outputs (Journel and Huijbregts, 1978; Goovaerts, 1997). In the machine learning literature, related models include multi-task GPs (Bonilla et al., 2008), semiparametric latent factor models (Teh et al., 2005), and multi-task kernel learning (Titsias and Lázaro-Gredilla, 2011). These can be viewed as variations of the LMC with different parameterizations and constraints. Convolution processes (CPs) have also been adapted to model multiple correlated outputs through the convolution of smooth kernels and latent processes (Álvarez and Lawrence, 2011)
. This approach usually has fewer hyperparameters and more efficient computation as compared to LMC, but only squared exponential (SE) kernels have been shown to be computationally tractable. Applying a multi-task GP (MTGP) framework(Bonilla et al., 2008) to clinical time series analysis has also been considered in two studies (Ghassemi et al., 2015b; Dürichen et al., 2015); both studies considered one patient as one task and used the remaining patients as reference training data. Other work adapted the LMC framework with one SE kernel to model three sparsely sampled clinical covariates (intracranial pressure, mean arterial blood pressure, and Pressure-Reactivity Index) jointly (Ghassemi et al., 2015b). The MTGP was shown to outperform a single-task GP (STGP) in prediction error. Both MTGP and CP have also been used with an SE kernel to model three densely sampled vital signs (respiration rate, systolic blood pressure, and heart rate); both methods showed improvements as compared to a single-task GP (Dürichen et al., 2015).
Our work is distinct from previous research in several ways. First, we use the GP regression framework to model multiple irregularly sampled medical time series using a sparse structured multi-output kernel. In contrast to related work (Ghassemi et al., 2015b; Dürichen et al., 2015), our kernel uses a mixture of flexible spectral kernels (Wilson and Adams, 2013), allowing periodic behavior and both short-term and long-term dependencies within and across the clinical covariates over time. Second, we use the LMC framework to enable an interpretable quantification of cross-correlation and sparsity between covariates. Third, we model many more clinical covariates (24) compared with previous studies (at most three); in the online medical setting, efficient and scalable computation in this multi-view model is essential. To the best of our knowledge, this is the first work that uses a sparse and low-rank formulation of the shared covariance matrix across clinical covariates to estimate and regularize the relationships between covariates in order to learn about covariate relationships specific to patient subgroups and to prevent overfitting.
In our methodology, MedGP, we trained a GP model on each reference patient separately, and used these models to estimate the empirical population-level model using nonparametric density estimation. This approach avoids training procedures that iterate through all reference patients, which is computationally intractable for an online system (Ghassemi et al., 2015b; Dürichen et al., 2015). To speed up training, we optimized the implementation in C++ using multithreading. Finally, in order to personalize the model for a new patient, we update the empirical population-level model on-the-fly to estimate patient specific parameters as measurements from the new patient are observed.
In this section, we describe our method, MedGP, for estimating the underlying dynamic processes jointly across a large number of sparsely sampled clinical covariates. We first describe the design of the Gaussian process kernel for capturing the temporal correlations within and between covariates. Next, we introduce the sparsity-inducing prior to regularize the LMC weight matrix. We then describe estimation of the parameters in the prior and the kernel. Next, we describe how to learn a patient-specific kernel by first building a population-level model from reference patients and then performing online updating of the parameters when observations about a new patient accumulate. Finally, we describe methods to perform computationally tractable online inference in these models, concluding with a discussion of computational complexity.
3.1 Gaussian Processes (GPs)
Gaussian processes (GPs) are distributions over arbitrary functions. By definition, a Gaussian process is a collection of random variables, any finite collection of which have a joint Gaussian distribution. Alternatively, a GP can be described as a distribution on an arbitrary function, defined as
where is the mean function:
and is the covariance function or kernel:
Any finite number of function values jointly have a multivariate Gaussian distribution with mean vector and covariance matrix between any pair of observations, defined by the kernel function,
Properties of the function such as smoothness or periodicity are determined by the kernel function . One of the most commonly used kernels is the squared exponential (SE) kernel
which is parameterized by a length scale and a scale factor . The functions generated by a GP with an SE kernel are smooth because the kernel function is infinitely differentiable (Rasmussen and Williams, 2006). The value of the length scale determines the distribution of changes over the function value with respect to changes in the input , encouraging a specific smoothness. Due to its simplicity, SE is used in many applications; however, the properties of the functions that it captures are fairly limited. Periodic functions, for example, are not well modeled by an SE kernel, but instead captured by a periodic kernel
where is the period of the function. When modeling medical time series, the SE kernel or the periodic kernel are often used in combination to capture the unknown source-specific smoothness and periodicity of the trajectories of clinical covariates (Stegle et al., 2008; Dürichen et al., 2015).
3.2 Gaussian Process Regression with a Structured Multi-Output Kernel
Our first goal is to jointly model multiple clinical covariates—vital signs and lab tests—over time for each patient using GP regression. For the th patient, we denote the time series of the th covariate as a vector , representing the time points that the th covariate was observed, and the corresponding observation vector :
where indexes time, and is the total number of observations for the th covariate of the th patient.
To represent the time series data over all covariates, we define
where , . Let be a multi-output function over time for the th patient. We capture the relationship between time and clinical observations as a GP regression model:
where is the residual noise vector. Marginally at the th observation of the th covariate, the residual noise is modeled as
is the covariate-specific residual variance for each individual.
We assume that the function is drawn from a patient-specific Gaussian process with mean function and kernel :
We set (Rasmussen and Williams, 2006).
We designed the kernel to capture predictive and generalizable covariance structure across medical time series data. Assuming the covariates are correlated across time, we adapted the linear model of coregionalization (LMC) framework (Journel and Huijbregts, 1978; Goovaerts, 1997). We used a set of basis kernels to model covariates jointly. The kernel for the cross-covariance of any pair of covariate types is modeled by a weighted, structured linear mixture of the basis kernels. The full joint kernel is written as a block structured function
where scales the covariance (defined by the th basis kernel) between covariates and , and . We collapsed into a set of weight matrices , where each is a symmetric positive definite matrix
If the inputs are the same for all covariates, we can further simplify Eq. (14) with Kronecker product . That is, if and :
although in practice we do not often see this situation in medical time series data.
The properties of the time series observations, such as periodicity and short term dependencies, are captured in the basis kernels. For medical covariates, the properties of each patent’s time series observations may vary. As a trivial example, when a patient is under age 18, their pulse will be well correlated with their age, height, and weight; above age 18, the correlation among pulse, age, height, and weight is more variable within age than across ages. Furthermore, only a few vital signs, such as heart rate, blood pressure, and body temperature, are known to be periodic with a 24-hour period (i.e., a circadian rhythm), but whether there is a similar period for specific lab results, such as white blood cell count or pressure of carbon dioxide in the blood, is unclear (Widmaier et al., 2004).
To handle the heterogeneity of patterns within covariates and across patients, we selected the spectral mixture (SM) kernel as the basis kernel (Wilson and Adams, 2013). The SM kernel is a general form of a variety of stationary kernels, including the squared exponential (SE) kernel and the periodic kernel, and has also shown good performance in modeling processes generated from more complex kernels through a mixture of kernels approach (Wilson and Adams, 2013). The basis kernel is written as
In our work, the mixture weights for each basis kernel are encoded in .
To be used for GP regression, must be a valid Mercer kernel, i.e., the Gram matrix must be positive definite for all and . Since the matrix produced by each basis kernel is symmetric positive definite, we only need to ensure that every is positive definite to produce a Mercer kernel. To do this, we parameterized as
Here , . We let denote the number of non-zero columns in , or the rank for when .
For any two observations from the same patient of different covariates at different times, denoted as and , the prior covariance from the GP kernel is
We summarize the parameters and hyperparameters of our SM-LMC kernel in Table 2.
|squared exponential part of th basis kernel|
|periodicity of th basis kernel|
|weights of for th basis kernel|
|intra-covariate weights of the th covariate for th basis kernel|
3.3 Sparsity-Inducing Priors on Weight Matrix
As the number of medical covariates included in the model increases, we need to increase the number of basis kernels and corresponding to allow greater representational flexibility. However, too many basis kernels may lead to overfitting and will become computationally intractable. To avoid this, we regularized the elements of each weight matrix by introducing structured sparsity-inducing priors on each matrix as follows.
We included two layers of sparsity-inducing priors for flexible, data-adaptive shrinkage behavior, modified from previous work (Polson and Scott, 2010; Gao et al., 2013). First, we put column-wise sparsity-inducing priors to regularize each column in
. This corresponds to regularizing the degree of freedom of the functions, or number of latent processes generated from each basis kernel in the LMC model(Álvarez et al., 2012). Second, we put sparsity-inducing priors on each matrix element in to produce element-wise sparsity. The effect of element-wise sparsity is to perform model selection on the number of basis kernels that each pair of covariates uses for covariance representation. Finally, we put sparsity-inducing priors on the elements of to shrink the covariance for observations from the same covariate.
In practice, we implemented each layer of the prior as a two-layer hierarchical gamma distribution. The generative model is written as
where each element has a Gaussian distribution. Parameters and control the column-specific shrinkage, while parameters and control the local shrinkage of each element in the matrix. For vector , we regularized each element with a local Laplace prior:
For our results, we set to recapitulate two layers of the horseshoe prior, using a statistically equivalent prior represented by a hierarchical gamma with four layers (Carvalho et al., 2010; Armagan et al., 2011; Gao et al., 2013; Zhao et al., 2016). Parameters , , , and were estimated during optimization. We set to regularize the diagonal terms . The hyperparameter controls the overall shrinkage profile of the hierarchical gamma prior (see Appendix A for more details). We chose over using cross-validation prediction error. Hyperparameter was chosen using grid search over the range using cross-validation prediction error as the objective.
3.4 Parameter Learning
To estimate the parameters for the regularized kernel, we optimized the posterior probability. We denote all parameters that were estimated directly asand hyperparameters in the sparsity-inducing prior as :
The posterior density of our model is then
The term is found by calculating the GP marginal likelihood given the values of (Rasmussen and Williams, 2006), which is
We thus estimated by solving the optimization problem:
See Eq. (31) in Appendix B for the derivation of .
Due to the conjugacy of the hierarchical gamma priors, we optimized parameters , , , directly using maximum a posteriori (MAP) estimates of their posterior distribution (or mean when the mode does not exist). Our optimization procedure then consists of two parts. In the first part, we used the update equations to estimate , , , and directly. In the second part, we estimated parameters , , , and using a scaled conjugate gradient method to find the local maximum, conditioned on , , , and . (Details can be found in Appendix B Eq. (32)–(35) and Eq. (36)–(40).) We iterated over the two steps until the change in reached the convergence criterion () or until the maximum number of iterations ().
3.5 Estimating the Population-Level Model and Online Updating
The GP with the structured kernel described above lets us model the patient-specific joint dynamics between covariates within the same patient. We now describe how we built a population-level empirical prior from a set of mixture kernels estimated from all training patients, and how we apply this empirical prior to a new patient.
To estimate the empirical priors across reference patients, we trained one GP kernel for each patient separately, and then we clustered and extracted the distribution of the basis kernels (defined by hyperparameters and ). The idea here is that, when we observe a set of estimated mixture kernels, we would like to understand the high-level properties of these mixture kernels shared across covariates and patients in the same patient group, and then estimate the distributions of these hyperparameters through observations of basis kernels belonging to this cluster. For instance, a circadian rhythm (24-hour periodicity) may be observed in some covariates for some patients, but the period across patients could vary within a range. Across the space of and the spectral kernels vary substantially (Figure 2a). For each basis kernel that was estimated, the characteristic period is and the length scale is (Wilson and Adams, 2013). There are different ways to define the features of a kernel. Here we used the temporal features of the learned kernels directly (Figure 2b). The temporal spacing of two adjacent points is one hour, and we use the kernel values within a window of length 72 hours. We then used a Gaussian mixture model (GMM) model to perform clustering on the kernels, and we chose the best number of kernel clusters () based on Bayesian information criterion (BIC). For the MedGP implementation, we adapted the open source scikit-learn package (Pedregosa et al., 2011). We used version 0.18.1, with ten random restarts, a maximum of 2,000 iterations, and allowing each mixture component to have its own covariance matrix.
For each identified kernel cluster, we estimated one set of parameters and for the basis kernel, and the weight coefficients—elements in matrices, computed using the matrices and vectors. We do this by building an empirical distribution using kernel density estimation (KDE) with a Gaussian kernel over the GP kernel hyperparameters assigned to that cluster. The bandwidth of the kernel density estimator was chosen based on Silverman’s “rule-of-thumb” (Silverman, 1986). We estimated each new parameter using density weighted means with the density from the univariate KDE as the weights. Note that if there were multiple kernels in a patient cluster, the estimated matrices were added based on the additive assumption of our kernel before aggregating to estimate the population-level kernel for that cluster. To allow online updating, we estimated the elements of the new empirical matrix and vector corresponding to each new
matrix using singular value decomposition (SVD). For the univariate GP regression, we did not use density weighted means because we found them to be unstable; instead we used a grid-based search to identify the hyperparameters with the highest kernel density estimates.
As the number of vital signs and lab measurements of a new patient accumulate, we update the hyperparameters to estimate a patient-specific kernel. Indeed, we update the kernel sequentially every time a new observation arrives. To do this in a computationally tractable way, we used the momentum method (Rumelhart et al., 1988) with a 72-hour window of previous observations to update the kernel hyperparameters when predicting the value of next observation. For all experiments, we chose the momentum as and the learning rate as . For elements in the matrices, we do not update the values if the elements were regularized to be zero so as to maintain the empirical sparsity structure.
3.6 Efficient Inference in MedGP
The main bottleneck of our method is in learning patient-specific kernel hyperparameters. Let denote the total number of samples of the th patient; the computational cost to compute the Gram matrix is , which increases linearly with the chosen number of basis kernels. To find the MAP estimates of the parameters, we need to invert and compute the determinant of the Gram matrix in Eq. (26). The computational complexity for the full matrix inversion is using Cholesky decomposition. When calculating the gradients for optimizing the hyperparameters, the cost is dominated by after the inverse Gram matrix is pre-computed, which is linear with the total number of the kernel hyperparameters. In practice, the complexity of each iteration is either or . That is, the patient with the most measurements is the main bottleneck for training. In our implementation, we mitigate the bottleneck using optimized linear algebra functions in Intel MKL library with multithreading and computing the gradients of the hyperparameters in parallel.
3.7 Medical Data Preprocessing
The HUP medical time series data consist of electronic health records (EHRs) from more than 260,000 patients admitted to a University of Pennsylvania Hospital. For each patient, the data include many heterogeneous clinical covariates, including ICD-9 codes, patient demography, length-of-stay, vital signs, and lab results. We jointly modeled the 24 covariates with the greatest number of observations across patients (Table 1). We selected three groups of discharged patients from these data: 1,365 septic patients, 952 patients with heart failure, and 4,723 patients with neoplasms. Each patient has at least one observation for each of the 24 covariates, and in total over four million observations were evaluated.
For each clinical covariate, we first removed obvious artifacts (e.g., values outside of the possible range in living humans). For the patients with neoplasms or heart failure, we used the full patient length-of-stay in training and testing. For septic patients, the disease progression varies substantially across patients, and the distribution of the covariates changes dramatically depending on the disease phase. To address this issue, we segmented the time series data into four disjoint partitions based on clinical status: no sepsis, pre-sepsis, sepsis, and recovery. To label each stage, we incorporated prior clinical domain knowledge. For instance, we identified sepsis stages using ICD-9 codes and positive blood culture results. Since our model assumes stationarity, to better estimate the temporal correlation across covariates, we chose the recovery stage before the patients’ discharge to test our method, since this is a relatively stable stage. We used the bed unit information to identify if the patient is in a stable state. That is, when a patient is transferred to step-down bed, we labeled the time series after the transfer as recovery. The median length-of-stay after pre-processing is 140 hours for the sepsis group, 285 hours for the heart failure group, and 197 hours for the neoplasms group.
We applied similar preprocessing procedure to the MIMIC-III data. We selected patients with a heart failure diagnosis that eventually had a routine discharge. We removed artifacts such as out of bounds values for each covariate, and applied the criteria to each patient that at least five measurements were taken for all 24 selected covariates. We extracted 1,004 heart failure under these criteria and used 1,003 of them, excluding one patient with more than 50K measurements due to memory constraints.
3.8 Experimental Setup
We applied MedGP to the three selected groups of patients separately, and evaluated characteristics and performance of MedGP under two different experimental settings. In the first analysis, we evaluated the model’s ability to learn the covariance between a pair of highly correlated clinical covariates, and we measured the imputation performance in an online setting. In the second analysis, we follow the same online setting, but instead jointly model all 24 clinical covariates, including four vital signs and 20 lab covariates. In both settings, we evaluated our method using 10-fold cross-validation at the patient level. That is, for each fold we ran the kernel clustering step on the kernels from the training patients to estimate a set of population-level basis kernels and matrices. This set of kernels was then applied to the held-out patients to predict the value of each covariate using observations from all other covariates measured at the same time as, or earlier than, the test observation (i.e., no future information included). After each prediction, we updated the patient-specific kernel parameters using the new observations from the test patient.
We compared our method to several univariate methods that modeled each covariate separately: (i) a naive one-lag prediction procedure, which predicts an observation equal to the last observation available from the same patient; (ii) an independent GP with squared exponential (SE) or spectral mixture (SM) kernels fitting each covariate separately (we tested with
for SM); (iii) the multi-resolution Probability Subtyping Model (PSM) combining linear regression, B-splines, and independent GPs(Schulam et al., 2015). To estimate the spectral kernel parameters, for each patient we initialized random kernels by drawing uniformly from a length scale range (between 6 and 72 hours) and period range (between 24 and 72 hours). We computed the marginal likelihood of all random kernels for each patient, and then initialize optimization using the kernels with the highest marginal likelihood. The elements in the matrices are initialized randomly between and .
We compared results from MedGP to these various methods using two metrics: (i) mean absolute error (MAE) of the predicted observations with the true observations, and (ii) 95% coverage, the percentage of true observations that fell within the predictive 95% confidence region. We quantified and reported the improvements with respect to both metrics compared to all three baselines (naive prediction, univariate GP, and PSM). To test if the differences in prediction results from different approaches were statistically significant, we performed paired t-tests for the results of each covariate and compared the-values with a Bonferroni corrected threshold (dependent on the number of jointly modeled covariates in each experiment).
We note that the original PSM was designed to model scleroderma disease (Schulam et al., 2015). Thus, to make it applicable to our different patient groups, several adjustments were made. First, we omitted the population and environmental factors selected for their relevance to scleroderma. Second, we chose the knots of the B-spline basis by sampling every hour for vital signs and every 24 hours for lab results between zero and the longest length-of-stay for patients in each disease group. Third, to make PSM training feasible on the scale of our data set, we limited the maximum number of subtypes to ten for the sepsis and heart failure groups, and 20 for the neoplasms group.
4 Experimental Results
We analyzed the performance of the method, MedGP—multi-output GP with a sparse SM-LMC kernel and online updating—by applying it to time series data from the Hospital of the University of Pennsylvania (HUP) and the public MIMIC-III data set (Johnson et al., 2016). We ran two types of experiments—one on two correlated lab covariates and the other with 24 covariates jointly. The results were compared against the baseline methods for prediction accuracy and 95% coverage calibration.
4.1 Results of Two Lab Covariates
As a proof of principle, we jointly modeled two well correlated lab covariates, prothrombin time (PT) and international normalization ratio (INR) on three HUP subgroups. PT measures the time it takes for the plasma in the blood to clot, and is often ordered to check bleeding problems. INR is an international standard for PT to account for possible variations across different labs. For the same patient, the two covariates usually have similar trajectories over time (Figure 1).
We trained the kernels for one patient’s INR and PT time series data both with and without the structured sparse prior (Figure 3). Both and matrices estimated using the sparse prior have higher levels of sparsity versus those estimated without using the sparse prior. We observed that for both methods, one of the estimated basis kernels captures long-term (around one month) dependencies. However, with the sparse prior, the estimated weights associated with this long term kernel are rank one instead of rank two. This means the trajectories of the two covariates are similar enough to be explained by one instead of two functions, and thus fewer hyperparameters. Moreover, two basis kernels were found with zeros weights and (Figure 3b), suggesting that the prespecified number of basis kernels may be reduced. We also found that the off-diagonal elements in the matrices in both cases have nonzero values, suggesting a nonzero covariance between PT and INR observations. In particular, two basis kernels captured the covariance between PT and INR: one with a greater than one-month trend (Figure 3b, and ), and one with a 27-hour trend (Figure 3b, and ). Here, the sparse kernel has 18 non-zero hyperparameters, whereas there are 40 for the non-sparse kernel. We can compare the two fitted kernels using both log marginal likelihoods and model selection scores. The log marginal likelihoods of the two kernels are (SM-LMC) and (sparse SM-LMC), indicating a better fit for the SM-LMC model without sparsity. However, the Bayesian information criterion (BIC) values, which take into account the number of parameters in a model, were (SM-LMC) and (sparse SM-LMC), where values closer to zero reflect better models. Thus, using a sparse prior has the advantage of a more compact kernel representation.
We then ran our model on all three disease groups separately, and compared our method with the univariate baselines described in Section 3.8 under the scenario of online imputation of the same two well-correlated clinical covariates. For independent GPs, we used gradient descent to optimize the hyperparameters. For PSM, we performed grid search for the parameters of the B-spline and the independent GP kernel. For our method, we set and for the matrices for training. In the sepsis and heart failure groups, three nonzero basis kernel functions () were found for the model using the SM-LMC kernel, while only two nonzero basis kernel functions () were found using the sparse SM-LMC kernel; the number of nonzero hyperparameters were and respectively. In the neoplasms group, the number of nonzero basis kernels were the same as the pre-specified number (). With 10-fold cross-validation, we found that results using the SM-LMC kernel showed smaller imputation error than those using the baselines for both PT and INR (Figure 4). The mean absolute errors (MAEs) showed that the non-sparse SM-LMC kernels perform imputation the best among the related approaches. On the other hand, looking at the 95% coverage, results using non-sparse or sparse SM-LMC kernel were well calibrated with respect to the confidence region compared with independent GPs, although sometimes slightly worse than PSM. Note that in this experiment we used a p-value threshold to detect statistical significance, which reflects the Bonferroni correction. The results indicate that the sparse prior finds models with sparse structure while maintaining prediction performance in this two covariate case.
4.2 Results of a Joint Model Including 24 Vital Signs and Lab Covariates
In the second experimental setting, we jointly modeled 24 vital signs and clinical covariates () for all three disease groups (Table 1). We set the number of basis kernels and the number of nonzero columns in as in this experiment for the three HUP subsets. For the MIMIC-III heart failure subset, we set . More detailed results of the best setup as well as the results with different could be found in Appendix C and Appendix D.
4.2.1 Estimating Population-Level Kernels
We first visualized the population-level kernels estimated from the three patient groups of the HUP data (Figure 5–7) and the MIMIC-III patient subgroup (Figure 8). We observed shared patterns in the basis kernels and the weight matrices across all patient groups. Comparing the estimated population-level kernels, we found at least one long-term smoothing basis kernel with length scale longer than three days, and one 24- to 25-hour periodic basis kernel, which indicates the existence of circadian rhythms in specific covariates as expected. Furthermore, in the neoplasms group, which consists of more patients than the other two groups, we found additional short-term smoothing basis kernels and one 12- to 13-hour periodic basis kernel, which may correspond to known circasemidian rhythm of clinical covariates, such as body temperature. We also observed an 11-hour periodic kernel in the MIMIC-III subset.
In addition to the characteristics of the basis kernels, our model with the sparse prior also showed interpretable cross-covariate patterns (Figure 5b, Figure 6b, and Figure 7b and Figure 8b). Based on the matrices, we identified groups of well correlated covariates. For instance, lab covariates hematocrit (Hct), hemoglobin (Hgb), and red blood cell (RBC) count showed the highest levels of correlation. Since both Hct and Hgb are known to be proportional to the number of red blood cells, this positive correlation was encouraging (Widmaier et al., 2004). The pair of lab covariates studied in the previous section, INR and PT, also showed substantial positive correlation. We found that the four vital signs—respiration rate (RR), heart rate (HR), systolic blood pressure (SBP), and body temperature (Temp)—had substantial correlations with each other as well as weak correlations with some lab covariates. Another identifiable set of well-correlated covariates includes lab measurements of carbon dioxide (CO), calcium, chloride, potassium, and sodium. The three lab covariates related to the concentration of hemoglobin—mean cell hemoglobin (MCH), mean cell volume (MCV), and mean cell hemoglobin concentration (MCHC)—appeared to have substantial correlation (Figure 5). The correlations modeled in these covariance matrices are exploited for accurate prediction and imputation in the MedGP framework.
To learn more about the importance of each kernel type across all subsets, we visualized the percent coverage of each type of kernel clusters found in all the subsets we have worked on (Figure 9). The coverage of each kernel type is computed as the ratio of patients that have non-zero matrix corresponding to it. We found that the kernel clusters of long-term (length scale days) and short-term (length scale hours) smooth dependencies have the highest coverage across four subsets. In the MIMIC-III subset, the coverages of the short-term kernel, and the 12-hour and 24-hour periodic kernels are higher than that of in the HUP subsets. We think this is because the higher sampling frequency in the MIMIC-III subset enables more accurate estimation of the short-term and periodic dependencies.
4.2.2 Results for Online Imputation
Next, we used the trained kernels to perform online imputation for each patient subgroup, where the goal is to predict the next observation for each covariate given the observations at previous time points. Across these methods, we used the percentage of improvement in MAE over three types of baselines—naive prediction, univariate GP (with SE or SM kernel), and PSM—to compare results for each of the 24 clinical covariates; we visualized the results separately (Figure 10–12; Figure 14–17 in Appendix C; Figure 18–33 in Appendix D). We also show the results of variations of our method for comparison (with or without the proposed sparse prior; with or without online updating). We performed paired t-tests on predictions from MedGP and each baseline to quantify the improvements, and statistical significance was evaluated using Bonferroni-corrected .
Comparing results with the independent GP model—specifically, selecting the best results from the SE or SM kernel, we found that MedGP, and in particular sparse SM-LMC with online updating, outperformed the independent GP model on the online imputation task for most covariates across the four patient groups (Figure 10). In the HUP data, we found 18, 21, 22 covariates significantly improved by MedGP in the sepsis, heart failure, and neoplasms subgroups respectively. In the MIMIC-III subset, we found 19 covariates were improved. For all four groups, the number of covariates that were improved significantly by MedGP is greater than using SM-LMC kernels without the sparse prior. We found that the covariates that were well correlated in usually showed significant positive improvements over independent GPs; Hct, Hgb, and RBC are notable examples. Similar observations could be made for INR and PT, the pair of lab covariates studied previously (Figure 4). Across 24 covariates, the MAEs for INR and PT were slightly worse compared with only modeling these two covariates. However, we also observed that using the sparse prior with the SM-LMC kernel led to better performance as compared to not using the sparse prior, indicating that sparse regularization is helpful when jointly modeling heterogeneous covariates. Finally, there were some covariates for which MedGP did not improvement over univariate GPs in two or more disease groups, including red cell distribution width (RDW), white blood cell count (WBC) and platelets.
When the baseline method is the naive one-lag method, for all three disease groups, we found fewer covariates with significant improvements compared with improvements over univariate GPs (Figure 11). In particular, the covariates for which the naive method had an advantage were lab covariates that have piece-wise linear behavior, such as mean cell hemoglobin (MCH) and mean cell hemoglobin concentration (MCHC) (Figure 1). In the case of piece-wise linear behavior, our kernel does not improve the performance compared with the naive approach since the time series are neither smooth nor periodic. Moreover, we also found that the naive method performed better in respiration rate, PTT, platelet, RDW, and white blood cell (WBC) count. Overall, however, our method improved online prediction results for 18, 20, 20 of the 24 covariates in sepsis, heart failure, and neoplasms groups, respectively. In the MIMIC-III subset, we found 14 covariates were improved significantly over the naive method.
When the baseline method is PSM (Schulam et al., 2015), we found that our method outperformed PSM for most of the lab covariates, but PSM outperformed MedGP in imputation of vital signs and two lab covariates: glucose point-of-care (Glucose POC) and potassium (Figure 12). For vital signs and glucose level, PSM has an advantage because of a higher sampling rate in those covariates and the highly structured mean function in PSM in the HUP subsets. The sampling rates are usually every 4 hours for vital signs and every 8 hours for glucose, which is more frequent than other lab covariates. Since PSM uses a B-spline basis function to capture the empirical mean, it may tolerate non-stationarity better. However, in the MIMIC-III subset, we observed that our method improved in imputing glucose and three vital signs (RR, SBP, temperature) over PSM significantly. We think this reflects the higher sampling rate of the covariates that allows better estimation of the short-term temporal dependencies. Overall, MedGP significantly improved the imputation of 17, 20, 18 covariates in sepsis, heart failure, neoplasms subsets respectively in the HUP data set, and 16 covariates in the MIMIC-III subset when compared with PSM. We contrast the PSM approach of structuring the mean function with our approach of structuring the kernel function, which leads to different types of gains in this problem.
Next, we looked at the calibration of the 95% coverage (Figure 16 and Figure 17 in Appendix C; Figure 26–33 in Appendix D). We found that MedGP outperformed independent GPs in terms of calibration of the 95% confidence region for all covariates. For this evaluation, the values closer to 95% are better. We observed that the coverage using the non-sparse SM-LMC kernel was usually higher than the coverage using the sparse SM-LMC kernel in the three HUP subgroups, indicating that MedGP may slightly underestimate covariate-specific noise. In contrast, in the MIMIC-III subset, we observed that MedGP gave consistently more accurate 95% coverage than without regularization in most covariates. We also found that, in all patient subsets, online updating significantly improves the accuracy of the 95% coverage. Among all tested methods, PSM tended to overestimate the 95% confidence region. We think this is because PSM assumes that the input time series are aligned by patient status, and this alignment is not the case in our data. With unaligned data, PSM learned large marginal variance parameters due to high empirical variance of the observations across patients at the same elapsed time. In contrast, the estimation of marginal covariance parameters in MedGP is not affected by alignment because estimates are patient-specific. We also observed that for either MedGP or PSM, the coverage was lower for some covariates in the MIMIC-III subset than in HUP subsets, such as temperature, CO, and PTT. This potentially reflects greater non-stationarity in the MIMIC-III subset, whose records were from intensive care units (ICUs) instead of regular hospital beds.
Finally, we compared the prediction performance of MedGP compared with the version without patient-specific online updating. We observed that online updating significantly improves the imputation errors of at least 12 out of 24 covariates in sepsis, heart failure, neoplasms, and the MIMIC-III subset (Figure 14 and Figure 15 in Appendix C; Figure 18–25 in Appendix D). Similarly, evaluating the 95% coverage, all 24 covariates were improved by the online updating across the three diseases groups in HUP, and 18 covariates were improved in the MIMIC-III subset (Figure 16 and Figure 17 in Appendix C; Figure 26–33 in Appendix D). This improvement highlights the importance of updating the empirical priors with patient-specific observations for this problem.
4.3 Computational Efficiency and Scalability
In this section, we compare computational speed between different implementations of our method. For patients with only a few observations, an existing implementation using conventional GP inference is sufficient for computationally tractable online inference. However, since our data include a large number of patients with potentially thousands of observations each, we implemented an exact inference algorithm in C++ and optimized it through Intel MKL libraries and customized multithreading blocks. In the experimental setting of , , and , there are 1114 hyperparameters to estimate. We summarized the runtime under different implementations for one patient with 2,028 unique time points and 6,679 observations (Table 3); the tests were performed using a machine with Intel(R) Xeon(R) CPUs running at 2.40GHz. Using our optimized implementation, for patients with large number of observations (), we accelerated training by a factor of 10 to 25 on average as compared with the sequential approach. We also compared our implementation with the standard GPy (GPy, since 2012) implementation under different sample sizes and , and reached empirically at least three times speed up. We provide these results in Appendix E.
|Computing Gram matrix||11||2|
|Inverting Gram matrix||13||3|
|Total per iteration||2521||102|
The proposed framework can be parallelized at the patient level and is suitable for analysis when patient data are observed in a streaming form. For each reference patient, we distributed the optimized training process on a computing cluster to estimate the patient-specific hyperparameters in parallel. In addition, the population-level kernels could be updated sequentially; the computationally expensive GP training procedure does not need to be applied to patient data in bulk. That is, when we receive more data from new patients, we only need to update the kernel density estimators. Our framework provides better computational efficiency compared to models designed for smaller collections of observations (e.g., approximately two hundred observations for each patient) as in most previous work. Those approaches are computationally intractable when working on a set of rich patient observations of the magnitude of the HUP data due to large matrix inversions and summing marginal likelihoods across patients at each iteration.
In this paper, we propose a flexible and efficient framework for estimating the temporal dependencies across multiple sparse and irregularly sampled medical time series data. We developed a model with multi-output Gaussian process regression with a highly structured kernel. We fit this model using an optimized implementation of exact GP inference to three different disease groups in the HUP medical data set and the MIMIC-III ICU data set. We showed that our method, MedGP, improves performance for online prediction of 24 clinical covariates as compared with independent univariate GPs, a naive method of propagating the previous observation, and an earlier state-of-the-art approach, PSM (Schulam et al., 2015). We found that, for well-correlated covariates, our method improves online imputation performance substantially over the related methods in most tested covariates. The improvements over the naive one-lag prediction and univariate GPs were significant in both vital signs and lab covariates. We found that PSM was, in general, better at predicting vital signs with more densely sampled observations. However, our approach does not require patient time series alignment and shows better calibration of the 95% confidence region as compared to PSM.
There are several directions that will be explored using the MedGP framework motivated by the present results. The first direction is to allow time-varying covariances by specifically modeling non-stationarity. Some possible approaches to explore include incorporating state-space models or change point detection (Adams and MacKay, 2007; Saatçi et al., 2010), and extending those methods to work on multivariate scenarios. Another direction of interest is to consider latent subpopulation-level structured kernels through multivariate medical time series. We expect that our results could be further improved through incorporating hierarchical methods with proper features or metrics to represent the differences between patients within the same disease group and across disease groups more carefully. For instance, the original PSM used three levels of hierarchy based on the subpopulations of patients with scleroderma, including population level, subpopulation level, and individual level. Our model may benefit from such an approach, but more efficient inference procedures are needed to train on our large data set (Feinberg et al., 2017). We should point out that this is possible through, for instance, deriving corresponding stochastic variation inference (SVI) algorithm. For example, previous work develops the SVI algorithm for semiparametric latent factor model (SLFM) with (Nguyen and Bonilla, 2014), which could be generalized to apply to MedGP.
For future applications, we will use the framework to monitor the health status of patients in a hospital setting and identify those patients at high risk for acute diseases in order to assist with decision making in treatment plans. Specifically, MedGP can impute latent state in patients at any time point, including confidence region around those estimates; this latent state can be used for a number of downstream analyses which require complete knowledge of patient state at specific time points. For instance, the changes of dynamics and temporal correlations between two vital signs have been found to be useful for disease detection given high-frequency regularly sampled time series (Nemati et al., 2012; Lehman et al., 2015). We demonstrated that MedGP accurately estimates the temporal correlations in the presence of sparse, unaligned time series data for up to 24 covariates, and we would expect to further associate the cross-covariate dynamics to more complicated diseases, such as septic shock (Henry et al., 2015), where the interactions of multiple covariates are jointly taken into consideration for diagnosis.
Appendix A. Details of the Hierarchical Gamma Prior
In this appendix, we provide more background and visualization for the hierarchical gamma prior we used for regularization. For the convenience, we use to denote the gamma function, and to represent a gamma distribution with shape parameter and rate parameter .
Following Proposition 1 in Armagan et al. (2011), for a random variable
drawn from a normal distribution with two-layered gamma priors on variance
is equivalent to the hierararchy
whereis given as
In Figure 13, we visualized the density of in Equation (29) for , and under different values of (Armagan et al., 2011). In this case, the prior distribution of is equivalent to a horseshoe prior, and can be interpreted as the shrinkage coefficient (Carvalho et al., 2010).
Appendix B. Details of Gradient Computation and Update Equations
In this appendix, the equations for the objective function during optimization, update equations for the parameters in the sparsity inducing prior and the gradients for the hyperparameters of the GP kernel are listed as reference.
The objective function to optimize for training one patient, , is
For update equations, we quoted from Zhao et al. (2016):
For partial gradients used for optimization:
Appendix C. Detailed Results of Imputation Error and 95% Coverage
We organized the detailed results of online imputation on all 24 covariates under the best number of basis kernel ( for HUP subsets and for the MIMIC-III subset) in Figure 14 to Figure 17. For Figure 14 and Figure 15, the mean absolute errors (MAEs) for each covariate is shown (in the original unit of measure). In Figure 16 and Figure 17, we showed the percentage for the prediction lied within the 95% confidence region (i.e. 95% coverage). We put markers in the figures to indicate the best among all methods, and the comparison of MedGP (sparse SM-LMC with online updating) against other methods. The statistical significance were tested using paired t-tests on patient-level results.
Appendix D. Results under Different Number of Basis Kernels
In this appendix, we showed more detailed results of the experiments using different number of basis kernels. We ran experiments with for on all four subsets. The results include all three subgroups in the HUP data set and the MIMIC-III heart failure subset. We visualized the results in Figure 18–33. We noticed that for most of the covariates, the imputation performance (both MAE and 95% coverage) improves as the number of increases. We also observed that the best number of varies across covariates under different metrics. For instance, for lab covariates INR and PT, we observed that setting or reduces MAE compared with , but the coverage still improves after . Allowing more numbers of basis kernels increases the flexibility for customization, but also increases complexity and thus the risk of overfitting for some covariates or patients. Overall for HUP subsets and