Adapting machine learning techniques to censored time-to-event health record data: a general-purpose approach using inverse probability of censoring weighting

1.1. Potential of predictive models trained using electronic health data

As an example, consider risk prediction models for cardiovascular disease and related outcomes (e.g., heart attack, stroke). Recent systematic reviews have described over 100 risk models produced between 1999 and 2009 [1, 2], including the well-known Framingham Risk Score [3], Reynolds Risk Score [4, 5], and the recent American Heart Association/American College of Cardiology pooled cohort equations [6]. Most risk prediction models have been estimated using data from homogenous and carefully selected epidemiological cohorts. These models often perform poorly when applied to diverse, present-day populations [7].

The increasing availability of electronic health data (EHD) and other sources of big biomedical data represents a key opportunity to improve risk prediction models. EHD, which consist of electronic medical records (EMRs), insurance claims data, and mortality data obtained from governmental vital records, are increasingly available within the context of large healthcare systems and capture the characteristics of heterogeneous populations receiving care in a current clinical setting. EHD databases typically include data on hundreds of thousands to millions of patients with information on millions of procedures, diagnoses, and laboratory measurement. The scale and complexity of EHD provide an excellent opportunity to develop more accurate risk models using modern machine learning techniques [8, 9, 10, 11, 12].

1.2. Challenge of right-censored data

EHD and other sources of big data in health care are not collected to answer a specific research question. In many datasets derived from EHD, the length of time that we are able to collect information on a particular subject is highly variable among subjects. Therefore, a large fraction of subjects do not have enough follow-up data available to ascertain whether or not they experienced the health outcome or adverse event of interest within a given time period (which we henceforth refer to as the event status). In the language of statistical survival analysis, the time of the adverse event (referred to as the event time) is said to be right-censored if the follow-up ends on a subject prior to her/him experiencing an event [13]. Unfortunately, fully supervised machine learning and classification methods typically assume that the event status is known for all subjects while in our setting the event status is undetermined for subjects whose event time is censored and who are not followed for the full time period over which one wants to make predictions (e.g., 5 years).

1.3. Existing techniques for right-censored data

To handle event statuses that are unknown due to right censoring, previous work has either proposed using preprocessing steps to “fill-in” or exclude observations with unknown event statuses or adapting specific machine learning tools to censored, time-to-event data.

In the later category, several authors including Hothorn et al. [14], Ishwaran et al. [15], and Ibrahim et al. [16] describe versions of classification trees and random forests to estimate the survival distribution. Lucas et al. [17] and Bandyopadhyay et al. [18] discuss the application of Bayesian networks to right-censored data. A few authors have considered applying neural networks to survival data but typically assume that the possible censoring and event times are few in number [19, 20]. Additionally, several have considered adapting support vector machines to censored outcomes by altering the loss function to account for censoring [21, 22]. These approaches are all based on modifying specific machine learning techniques to handle censoring, which limits the generalizability of the approach used to handle right-censored outcomes. For example, to adapt decision trees and random forests to right-censored outcomes, the authors cited above modify the splitting criterion to accommodate censoring. Instead of splitting the data to minimize the node impurity, they choose the split which maximizes the log-rank statistic, a statistic to compare the difference in the survival curves between two groups (in this case between two child nodes). However, the recursive partitioning approach of decision trees is fundamentally different than, e.g., the approach of Bayesian networks or generalized additive models, so the idea of altering the splitting criterion in a tree to accommodate censoring does not apply to these other approaches.

Alternatively, several general-purpose ad hoc techniques have been proposed for handing observations with unknown event status including 1) discarding those observations [23, 24], 2) treating them as non-events [25, 26], or 3) repeating those observations twice in the dataset, one as experiencing the event and one event-free. Each of these observations is assigned a weight based on the marginal probability of experiencing an event between the censoring time and the time the event status will be assessed [27]. These simple approaches are known to induce bias in the estimation of risk (i.e., class probabilities) because, as we discuss later, discarding observations with unknown event status or treating them as non-events necessarily over- and under-estimates the risk [25, 26]. Even the third approach, although more sophisticated, produces poorly calibrated risk estimates because these weights attenuate the relationship between the features and outcome.

In summary, previous approaches for handling right-censored time-to-event data are specific to a single machine learning technique or are generally applicable but produce poorly calibrated risk estimates.

1.4. Our approach

The goal of this paper is to propose a general-purpose technique for mining right-censored time-to-event data which has improved calibration compared to the ad hoc techniques previously proposed. Specifically, we introduce a simple, pre-processing step which re-weights the data using inverse probability of censoring (IPC) weights. The IPC-weighted data can then be analyzed using any machine learning technique which can incorporate observation weights. Briefly, subjects with unknown event status are given zero weight; subjects with a known event status are given weights to account for subjects who would have had the same event time but were censored. Subjects with larger event times are assigned higher weights to account for the fact that they are more likely to be censored prior to experiencing the event of interest. This heuristic explanation is made mathematically precise later.

There has been some prior work which used inverse probability of censoring weighting (IPCW) in machine learning methods. For example, Bandyopadhyay et al. [18] discuss how to use IPCW specifically in the context of estimating Bayesian networks with right-censored outcomes. However, to the best of our knowledge, this paper is the first to propose the use of IPCW as a general-purpose technique that may be used in conjunction with many machine learning methods. IPCW properly accounts for censoring and can be easily integrated into many existing machine learning techniques for class probability estimation, allowing new machine learning tools for censored data to be created.

2.1. Notation and terminology

Let E be the indicator that a health outcome or adverse event occurs within τ years of a pre-defined timepoint. Throughout the paper, we refer to E as the τ -year event status. In our setting, for example, we are interested in whether or not a major CV adverse event (hereafter referred to as a CV event) occurs within 5 years of an “index” clinic visit where risk factor data are available. Binary classification methods typically assume that E is fully observed for all patients, but this is unlikely to be true when using current EHD. When a patient leaves the health system or the study ends before τ years of follow-up and before the subject experiences the event of interest, the subject’s event status at τ is unknown. For example, assume a subject is enrolled in the health system for only 3 years and does not experience a CV event. We do not know if this subject would have experienced a CV event within 5 years as we do not get to observe any health outcomes between 3 and 5 years. For this subject E would be unknown.

For individual i define T_i as the time between baseline and the event of interest, and define C_i as the time between baseline and when the patient is lost to follow-up (e.g., in our context, disenrolls from the health plan or reaches the end of the data capture period without experiencing an event). We observe V_i = min(T_i, C_i) and δ_i = 𝕀(T_i < C_i), the indicator for whether or not an event occurred during the follow-up period. If δ_i = 0, the i^th subject’s event time is said to be right-censored. The value of E_i is unknown if subject i is censored prior to τ, or equivalently if min(T_i, τ) ≥ C_i. Continuing the example described above we have V_i = 3, δ_i = 0, and since this subject is censored prior to 5 years, E_i is unknown. We will denote the set of features available on individual i by X_i; it is assumed that these features are fully observed at the beginning of the follow-up period and, hence, are not subject to censoring and do not vary over time. The target of prediction is π(X_i) = P(E_i = 1|X_i) ≡ P (T_i ≤ τ |X_i), and predictions are denoted by $\widehat{\pi }\left({\mathbf{X}}_i\right)$.

2.2. The IPCW method

We propose to use an inverse probability of censoring weighting (IPCW) approach for censored event times which is well-established in the statistical literature but has not been broadly applied for machine learning. Intuitively, excluding subjects for whom E is unknown leads to poor risk prediction because subjects with small event times are less likely to be censored than those with event times beyond τ. Therefore, we oversample subjects with E = 1 if we exclude patients for whom E is unknown. In IPCW, only those subjects for whom E is known contribute directly to the analysis, but they are reweighted to accurately “represent” the subjects with unknown E.

The advantage of IPCW to account for censoring is that it is a general-purpose approach that may be applied to any machine learning method. The analyst can then apply several different machine learning methods for risk prediction and select the optimal one based on censoring-adjusted criteria discussed in Section 4. The general-purpose IPCW method proceeds as follows:

1. Using the training data, estimate the function G(t) = P(C_i > t), the probability that the censoring time is greater than t, using the Kaplan-Meier estimator of the survival distribution (i.e., 1 minus the cumulative distribution function) of the censoring times [13],

   $$\widehat{G}\left(t\right)=\prod _{j:V_j<t}\left(\frac{n_j-d_j^{\ast }}{n_j}\right)$$

   where $d_j^{\ast }$ is the number of subjects who were censored at time V_j, and n_j is the number of subjects “at risk” for censoring (i.e., not yet censored or experienced the event) at time V_j. We note that, for IPCW, Kaplan-Meier is applied to estimate the distribution of censoring times, whereas it is much more commonly used to estimate the distribution of event times.
2. For each patient i in the training set, define an inverse probability of censoring weight,

   $${\omega }_i=\{\begin{array}{cc}\frac{1}{\widehat{G}\left\{\text{min}(V_i,\tau )\right\}}\hfill & \text{if}\text{min}(T_i,\tau )\le C_i\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\phantom{\}}$$

   (1)

   Patients whose event status is unknown at τ (i.e., are censored prior to τ and therefore have C_i ≤ min(T_i, τ)) are assigned weight ω_i = 0, and hence are excluded from the analysis. The remaining patients are assigned weights inversely proportional to the estimated probability of being censored after their observed follow-up time.

3. Apply an existing prediction method to a weighted version of the training set where each member i of the training set is weighted by a factor of ω_i. In other words, if ω_i = 3, it is as if the observation appeared three times in the data set.

Step 3 is left purposefully vague, as the specific manner in which IPC weights are incorporated will vary across machine learning techniques. Off-the-shelf implementations of some techniques allow for the direct specification of observation weights, in which case little additional work is needed to get risk estimates. More generally, most machine learning techniques involve estimation (typically using maximum likelihood estimators) and assessing model fit/purity and incorporating weights in both steps is straightforward. Section 3 illustrates how a variety of machine learning algorithms can be adapted to handle weighted observations and hence be “adapted” for censoring using IPCW.

2.3. Intuition of IPCW and a toy example

We briefly argue heuristically why inverse probability of censoring weighting appropriately handles censoring and leads to accurate risk prediction across a variety of machine learning techniques. Suppose we estimate that 1/3 of subjects have censoring times greater than 2.5 years (i.e., $\widehat{G}$(2.5) = 1/3), and that the i^th subject is observed in our study to experience an event at t = 2.5 years (i.e., δ_i = 1 and V_i = 2.5). For this subject, the event status is known (E = 1) and her/his IPC weight is ω_i = 3. This subject is weighted by a factor of 3 because she/he can be thought of as representing 3 individuals: 2 similar or “shadow” subjects censored prior to their event time at t = 2.5 for whom E is unknown, plus themselves (recall that on average 2/3 of subjects in this example with event times equal to 2.5 are censored prior to experiencing the event). Thus, subjects with known event status E and a longer time-to-event receive larger weights as they represent a greater number of “shadow” subjects whose event status is unknown due to censoring. IPCW is conceptually equivalent to creating a new dataset where each subject is replicated ω_i times. However, creating such an expanded dataset is often not advisable, both for reasons of practicality (memory/storage limitations) and mathematical precision (ω_i may not be an integer or simple fraction). A full justification of the use of IPCW can be found in Tsiatis [29]. We provide the main mathematical details in the Supplementary Material.

As an example, consider the following toy dataset given in Table S1 (see Supplementary Materials) with only 50 observations and a single binary covariate. Suppose that we wish to estimate the probability of having an adverse event within 5 years within each level of the covariate (i.e., τ = 5). If we knew E_i for all subjects, we would just take the average of E_i within each level of the covariate. However, because some subjects do not have 5 years of data and did not experience an adverse event during their follow-up, E_i is unknown for them (which is indicated by a question mark in the table).

Instead, to implement IPCW, we estimate the censoring distribution G(t) using a Kaplan-Meier estimator and compute $\widehat{G}${min(V_i, τ)}. The weight, ω_i for each subject is given by Equation 1. Now to estimate probability of having an adverse event within 5 years within each level of the covariate we take a weighted average of E_i within each level of the covariate; i.e. $\widehat{P}(E=1\mid \mathbf{X}=1)=\frac{{\sum }_{i:X_i=1}{E}_i{\omega }_i}{{\sum }_{i:X_i=1}{\omega }_i}=0.58$ and similarly $\widehat{P}(E=1\mid \mathbf{X}=0)=\frac{{\sum }_{i:X_i=0}{E}_i{\omega }_i}{{\sum }_{i:X_i=0}{\omega }_i}=0.53$. Subjects for whom E_i is unknown have weight equal to 0 so E_iω_i = 0. Of course, many machine learning methods are more sophisticated but the basic idea presented here is still applicable.

3.1. Logistic regression and generalized additive logistic regression

Logistic regression is a simple and popular technique for modeling binary outcome data. The goal is to find a linear combination of features to approximate the log-odds, i.e.,

where π(x) = P(E = 1|X = x) for the vector of features X. In risk prediction, the model often only includes the main effects of each risk factor, i.e., the value of the risk factor itself. Given features X and a corresponding vector of event indicators E, the logistic regression log-likelihood takes the form

where n is the number of observations in the training set and β = (β₀, …, β_p)^T. This log-likelihood can be maximized using a variety of techniques to find the maximum likelihood estimates.

In the case that E_i is unknown due to right censoring, IPC-weighted logistic regression simply maximizes the weighted log-likelihood

where contribution of the i^th subject to the likelihood in (2) is weighted by ω_i given in (1).

Logistic regression with only the main effects of various risk factors is unlikely to produce a well-fitting model when the log odds of experiencing the event has a non-linear relationship with the features. Enlarging the feature set by considering a basis expansion of the continuous features may improve prediction. If z_j is the basis expansion of the j^th feature and β_j is vector of the same dimension as z_j, then the generalized additive logistic model (GAM) assumes that

Restricted cubic smoothing splines, B-splines, or thin-plate regression splines are frequently used as the basis expansion in practice. Since expanding the feature space involves estimating many more parameters, it is common to penalize the roughness of the linear predictor ${\sum }_{j=1}^p{\beta }_j^T{\mathbf{z}}_j$, and estimate the regression parameters by maximizing the resulting penalized log-likelihood

where $\beta ={({\beta }_0,{\beta }_1^T,\dots ,{\beta }_p^T)}^T$, S_j, j = 1, …, p, are appropriately chosen smoothing matrices, and λ_j, j = 1, …, p are tuning parameters which control the degree of penalization/smoothness. λ_j are typically selected to minimize the unbiased risk estimator (UBRE), which in the case of logistic regression is proportional to the Akaike Information Criterion given by AIC = 2k − 2ℓ, where ℓ is the (log-)likelihood given in (2) and k is the total number of free parameters.

Similarly, the IPC-weighted generalized additive logistic regression maximizes the following weighted version of the penalized log-likelihood given in (4):

The scores used to select the tuning parameters in the generalized additive model are also easily modified using IPC-weights to account for right censoring. In particular, the weighted AIC becomes AIC^ω = 2k−2ℓ^ω, with ℓ^ω given in (3).

3.2. Bayesian networks

Bayesian networks have been used extensively in biomedical applications to aid in understanding of disease prognosis and clinical prediction [30, 31] and guide the selection of the appropriate treatment [32, 33] in clinical decision support systems. Lucas et al. [17] provide a comprehensive review of Bayesian networks in medical applications.

The key to Bayesian network techniques is that using Bayes theorem one can rewrite π(x) as

so that focus is now shifted to estimation of the conditional density/probability P_X|E(x|e) and the probability P_E(e) for e = 0, 1. When E is observed on all subjects (i.e., there is no censoring), the maximum likelihood estimate of P_E(e) is given by the sample mean of the event indicators. To simplify the task of modeling P_X|E, one can represent the joint distributions of X|E using a directed acyclic graph (DAG), i.e., a Bayesian network. One advantage of the Bayesian network approach is that clinical knowledge and data can be combined to suggest and refine DAG structures. The DAG encodes conditional independence relationships between variables, allowing the joint distribution to be decomposed into a product of individual terms conditioned on their parent variables [34]:

where Pa(X_j) are the parents of X_j. Several approaches have been proposed to modeling the terms P_Xj|Pa(Xj),E{x_j |Pa(x_j), e}. In many applications, continuous covariates are discretized to allow learning of the joint density of P_X|E(x|e) nonparametrically. In the application considered in this paper, all parent nodes are discrete (or have been discretized) which simplifies the modeling considerably. If the j^th feature X_j is discrete, then P_Xj|Pa(Xj),E{x_j|Pa(x_j), e} is estimated by computing the proportion of observations in each unique state of X_j separately for each level of Pa(X_j) and each level of E via

which is the non-parametric maximum likelihood estimator.

To fit the Bayesian network using IPCW, we make the following modifications as discussed in Bandyopadhyay et al. [18]. We can obtain the IPCW maximum likelihood estimator of P_E(e) using the weighted mean ${\widehat{P}}_E\left(e\right)=\frac{1}{n}{\sum }_{i=1}^n{\omega }_i\mathbb{𝕀}(E_i=e)$. We note that this is equivalent to the Kaplan-Meier estimator of P_E (e). Similarly, we can then obtain an IPCW maximum likelihood estimator of the distribution for the discrete variables X_j separately for each level of Pa(X_j) and each level of E so that (6) becomes:

A number of parametric and semi-parametric approaches to modeling the covariate distributions of continuous features are possible and have been described elsewhere [35, 36], and we discuss in the Supplementary Material how to adapt these approaches for censored outcomes.

We note that tuning a Bayesian network for optimal performance may involve determining the network structure and/or controlling model complexity for a given structure. In the Bayesian network implementation for our data application, we consider only a single network structure which is informed by discussions with our clinical colleagues (see Figure 3); however, a set of feasible structures could easily be compared on a test set or via cross-validation using the calibration and reclassification metrics described in Section 4.

3.3. Decision trees

Recursive partitioning is a powerful and flexible way to build predictive models for both discrete and continuous outcomes, and decision tree algorithms are widely applied in biomedicine [see 37, 38, and references therein]. Decision trees aim to partition training data into subgroups with homogeneous outcomes, with subgroups defined by a set of binary splits of the features. The prediction for a given test instance is made by identifying the partition or node it belongs to, then computing a summary statistic (e.g., the sample average) for training instances in that partition or node.

Many techniques have been proposed to grow decision trees, mostly differing in the criteria used to decide how/if to split a node and to prevent overfitting. One popular technique, CART [39], uses the decrease in Gini impurity to determine which feature and at what level to split a node. The change in Gini impurity for a possible splitting rule is given by

where $\widehat{\pi }\left(s\right)$, $\widehat{\pi }\left(s_l\right)$ and $\widehat{\pi }\left(s_r\right)$ are respectively the sample proportion of outcomes in a node s, the node’s left-hand children s_l, and the node’s right-hand children s_r for the particular splitting rule; N_sl and N_sr are the number of instances in each child; and N_s = N_sl + N_sr. Alternatively, the C4.5 and C5.0 decision trees [40] use the information gain metric instead of the Gini impurity to make decisions on how to split each node.

In the unweighted case, the sample proportions $\widehat{\pi }\left(s\right)$ for node s are computed as the average of the event indicators E for subjects in node s. For a test instance with features x falling in terminal tree node s_T, we can estimate the risk π(x) as the proportion of training instances in that node with E_i = 1.

It is straightforward to extend decision trees to incorporate IPCW: individual cases in the training set are assigned weights ω_i as described above, and the ω_i are used as “case weights” in the decision tree algorithm. With IPCW, we calculate a weighted decrease in Gini impurity,

where

and

A similar approach could be applied to estimate a weighted version of the information gain metric.

Once the structure of the tree has been determined, the predicted risk of a test instance with features x falling in terminal node s_T is estimated using the weighted mean:

where ω_i is the IPC weight given in (1) and $N_{s_T}^{\omega }={\sum }_{i=1}^n{\omega }_i\mathbb{𝕀}(S_i=s_T)$.

Because of their flexibility, classification trees often overfit training data. Many overfitting avoidance techniques have been proposed, with most involving a tuning parameter which restricts the complexity of the tree. One strategy consists of setting a lower limit m on the number of individuals assigned to a terminal node; in our notation above, the node S would not be split according to a given rule unless min(N_sl, N_sr) ≥ m. This strategy is easily generalized to the case with censoring by requiring that min($N_{s_l}^{\omega }$, $N_{s_r}^{\omega }$) ≥ m. However we note that N_sl ≈ $N_{s_l}^{\omega }$ and N_sr ≈ $N_{s_r}^{\omega }$ as the expected value of ω_i is one, so in practice setting a lower bound for min(N_sl, N_sr) is usually sufficient. Another approach only pursues splits where the change in Gini impurity exceeds a certain threshold θ, e.g., ΔI_G(s) ≥ θ. Substituting Δ$I_G^{\omega }$(s) for ΔI_G(s) allows the same rule to be used in the censored data setting. Final tuning parameter values may be chosen by cross-validation, where the cross-validated criterion to optimize could involve a measure of calibration, discrimination, or a combination of both as discussed in Section 4.

3.4. k-nearest neighbors

The k-nearest neighbors classifier is widely used in biomedical applications and provides a flexible, powerful, and intuitive method for risk prediction [as examples, see 41, 42, 43, and references therein]. Define d(X_i, X_j) to be a distance metric between two vectors of features X_i and X_j. To estimate the event probability for an instance in the test set with features X = x, define R_i(x) to be the rank of the distance between x and X_i, i.e., d(x, X_i), among all n observations in the training data set. When the event status is known on all subjects in the training dataset, in the most straightforward application of k-nearest neighbors, π(x) is simply the proportion of the training instances experiencing the event among those with R_i ≤ k, the non-parametric maximum likelihood estimator.

The key choice for implementing a k-nearest neighbor classifier is to select an appropriate distance metric and the number of neighbors to consider. The number of neighbors may be treated as a tuning parameter and chosen using cross-validated estimates of some appropriate criterion discussed in Section 4.

To adapt a k-nearest neighbor classifier to the situation when E_i may not be known for all subjects, we note that the distance between the vector of features is not affected by censoring. Therefore, to predict the probability of an event for an instance in the test set, we can identify the k closest neighbors in the training set just as we did before. However, we now replace the simple average with a weighted one:

In this extension of the k-nearest neighbor classifier, we choose the k neighbors regardless of the value of the IPC weights, ω_i, for those neighbors in the training set. Therefore, the total sum of the weights for the k neighbors ${\sum }_{i=1}^n{\omega }_i\mathbb{𝕀}\{R_i\left(\mathbf{x}\right)\le k\}$ may be different and the number of training instances with non-zero weight among those k neighbors may vary depending on the features of the test instance. But since the expected value of the IPC weight is equal to one, independent of the features X_i of the training instance, ${\sum }_{i=1}^n{\omega }_i\mathbb{𝕀}\{R_i\left(\mathbf{x}\right)\le k\}$ should be approximately equal across different values of x, and we do not have to worry about adjusting the number of neighbors across the feature space.

3.5. Other machine learning techniques

We note that the manner in which the IPC-weights are incorporated in the estimation procedures (e.g., maximizing a weighted objective function) and the calculation of measures of model fit (e.g., weighted Gini impurity) in these four illustrative examples may be easily applied to other machine learning techniques. For instance, given an implementation of IPCW decision trees, constructing an IPCW random forest is straightforward. Furthermore, IPCW can also easily be incorporated into other, distinct machine learning methods such as support vector machines (SVMs) [44] and multivariate adaptive regression splines (MARS) [45] using the framework and implementation developed in Sections 3.1 – 3.4.

For example, SVMs seek to find the hyperplane in the feature space (possibly including basis transformations) that separate those that experience the event and those that do not by the largest distance while bounding the proportional amount by which some predictions are on the “wrong-side” of the margin. When E_i is known on all subjects, the SVM solution can be found from regularization function estimation [46], i.e., given by the solution to

where [x]₊ is 0 if x < 0 and equal to x if x > 0, λ is a tuning parameter that could be selected using cross-validation, z_i is the feature vector for the i^th subject including any basis transformations, and E_i is coded as +1 and −1 for subjects experiencing and not experiencing the event, respectively. The SVM classifier is given by sign$({\widehat{\beta }}_0+z^T\widehat{\beta })$.

When E_i is possibly unknown due to right censoring, we can easily solve an IPC-weighted version of Equation 11 just as we used an IPC-weighted objective function (i.e., the IPC-weighted log likelihood) in Section 3.1. That is, to obtain ${\widehat{\beta }}_0$ and $\widehat{\beta }$ using IPCW we could solve

Similarly, the framework discussed in Section 3.1 for generalized additive models could be used to adapt MARS when the event status may be unknown due to right censoring. Specifically, MARS uses pairs of expansions in piecewise linear basis functions which take the form [x_ij − t]₊ and [t − x_ij]₊ where the possible knots t are the observed values of the j^th covariate across all subjects in the training set. When the outcome is binary, a linear regression model fit using least squares is used to select the basis functions for inclusion in the model. Typically, a forward stepwise feature selection procedure is used in which the basis function pair (or product of a basis function pair with another term already included in the model) is selected to minimize the residual squared error. Once the basis functions have been selected, a logistic regression model is fit with those covariates [47].

Again, when E_i is unknown due to right censoring, each step of the process can be altered to incorporate IPC-weighting. Specifically, the linear regression model can be fit using IPC-weighted least squares; the basis function pair is selected to minimize the IPC-weighted residual squared error, and then an IPC-weighted logistic regression model is fit with the selected covariates (similar to Equation 3).

4.1. Calibration

In standard risk prediction settings, calibration is commonly assessed by ranking the predicted risks $\widehat{\pi }\left({\mathbf{x}}_i\right)$, partitioning the ranked predictions into bins B₁, B₂, …, B_m (e.g., by decile or clinically relevant cut points), and comparing the average predicted risk in each bin to an empirical estimate of the risk within that bin. When E_i is known for all subjects, the empirical risk estimate for bin B_k is simply given by ∑_i∈Bk E_i/|B_k|, where |B_k| is the number of instances in bin B_k. However, when the outcome E_i is unknown for some subjects within a bin, an alternative estimator of the empirical risk is needed. One option estimates the probability of experiencing an event prior to time τ within each bin using the Kaplan-Meier estimator, yielding a calibration statistic of the form:

where ${\stackrel{‒}{\pi }}_k$ is the average of predicted probabilities in bin k, ${\widehat{p}}_k^{KM}$ is the Kaplan-Meier estimate of experiencing an event before τ among test subjects in bin k, and $var\left\{{\widehat{p}}_k^{KM}\right\}$ is the sampling variance of the Kaplan-Meier estimator calculated. K is analogous to the χ² statistic with m − 2 degrees of freedom for assessing the calibration of logistic models [48].

4.2. Concordance index

The area under the ROC curve (AUC) is a widely used summary measure of predictive model performance. When the outcome is fully observed on all subjects, it is equivalent to the concordance index (C-index), the probability of correctly ordering the outcomes for a randomly chosen pair of subjects whose predicted risks are different. As described in Harrell [49], the C-index can be adapted for censoring by considering the concordance of survival outcomes versus predicted survival probability among pairs of subjects whose survival outcomes can be ordered, i.e., among pairs where both subjects are observed to experience an event, or one subject is observed to experience an event before the other subject is censored. The C-index adapted for censoring is given by

where 𝕀(·) is the indicator function.

5.1. Example dataset: Predicting cardiovascular risk using electronic health data

We illustrate the application of IPC-weighted risk prediction methods to the problem of predicting the risk of a cardiovascular event from electronic health data. The data come from a healthcare system in the Midwestern United States and were extracted from the HMO Research Network Virtual Data Warehouse (HMORN VDW) associated with that system [52]. The VDW stores data including insurance enrollment, demographics, pharmaceutical dispensing, utilization, vital signs, laboratory, census, and death records. This healthcare system includes both an insurance plan and a medical care network in an open system which is partially overlapping. That is, patients of the insurance plan may be served by either the internal medical care network and or by external healthcare providers, and the medical care network serves patients within and outside of the insurance plan. Patient-members who do not visit any of the clinics and hospitals in-network do not have any medical information (e.g., blood pressure information) included in the electronic medical record (EMR) of this system. Furthermore, once the patient-member disenrolls from the insurance plan, the patient is right-censored as there is no longer any information on risk factors or outcomes (i.e., CV events) recorded in the EMR or insurance claims data.

This study and the use of these data were approved by the Institutional Review Boards of both the University of Minnesota and HealthPartners Institute for Education and Research.

5.2. Defining the study population

The study population was initially selected from those enrolled in the insurance plan between 1999 and 2011 with at least one year of continuous insurance enrollment and prescription drug coverage.

We included only patients with at least two medical encounters in the in-network ambulatory clinics (excluding urgent care) which had blood pressure recorded and were at least 30 days but at most 1.5 years apart. Patients who are only treated in the emergency room or urgent care clinics (i.e., settings where patients are unlikely to be counseled about their CV risk) were excluded in this analysis. As is typical in other CV risk prediction models [3], patients under the age of 40 and with pre-existing serious comorbidities other than diabetes (e.g., prior CV event, chronic kidney disease, etc.) were also excluded. After applying the above criteria, our final analysis dataset contained 87,363 individuals.

The available longitudinal data on each patient-member was divided into: (i) a baseline period, where the risk factors were ascertained, and (ii) a follow-up period, where we assessed whether a patient experienced a CV event (and, if so, when). The baseline period consisted of the time between the first blood pressure reading during the enrollment period and the date of the final blood pressure reading at most 1.5 years from the first measurement. Although some clinics in this health system adopted the electronic medical record as early as 1999, most clinics transitioned to the electronic system between 2001 and 2002, so the earliest blood pressure readings recorded in the medical record are typically between 2001 and 2002. The follow-up period for a patient begins at the end of the baseline period, referred to as the index date, and continues until either the patient experiences a CV event (defined below), the patient disenrolls from the insurance plan, or the data capture period ends (in 2011), whichever comes first. The distribution of the follow-up periods for the resulting analysis cohort is shown in Figure 1, which illustrates that a large proportion of subjects’ CV event times are censored prior to the end of follow-up. Figure 2 shows that, unless we consider a very short time horizon τ, the τ -year event status will be unknown for a substantial proportion of subjects in this cohort. In particular, for τ = 5 years, the proportion of subjects for whom the τ -year event status is known is only 47.8%.

5.3. Risk factor ascertainment

Risk factors used as features in the machine learning models included age, gender, systolic blood pressure (SBP), use of blood pressure medications, cholesterol markers (HDL and total cholesterol), body mass index (BMI), smoking status, and presence/absence of diabetes. These risk factors were chosen because they have been consistently used in prediction of adverse cardiovascular outcomes in the work cited in Section 1.1. Summary statistics and brief descriptions for the risk factors are given in Table 1. Missing risk factor values were filled in prior to model fitting using multiple imputation by chained equations [53].

Unless otherwise noted, SBP, HDL, total cholesterol, and BMI were all considered as continuous features. For systolic blood pressure, we took the average of all the blood pressure measurements during the baseline period excluding readings obtained during emergency department visits, urgent care visits, hospital admission, and during procedures (e.g., surgeries) because they may be influenced by acute conditions. Use of SBP medication during the baseline period was inferred from claims data. Diabetes was defined based on joint consideration of inpatient and outpatient ICD-9-CM diagnosis codes (ICD-9-CM codes 250.xx), use of glucose-lowering medications, and glucose-related laboratory values using a previously validated algorithm with estimated sensitivity of 0.91 and positive predictive value of 0.94 [54]. Body mass index was calculated as a function of patient’s height and weight which are recorded in the EMR. The height of an individual is the average height measured at any encounter (possibly outside of the baseline period). Because all subjects in the analysis dataset are adults, we expect height to remain relatively constant over the follow-up period. The weight is calculated as an average of all weight measurements taken during the baseline period. In the EMR in this health system, smoking status is categorized as never smoked, smoking, quit smoking, and passive (i.e., second-hand) smoking. In our analysis, a person is considered to have never smoked only if they consistently recorded “no smoking” throughout the baseline period. For the purpose of constructing the model we combine the “passive smoking” and “no smoking” categories. Finally, the most recent laboratory measurements before the end of the baseline period for HDL and total cholesterol were used.

5.4. Cardiovascular event definition

Cardiovascular events were defined as the first recorded stroke, myocardial infarction (MI), or other major CV event after the baseline period, prior to 5 years of follow-up. Major CV events were ascertained based on the date of primary hospital discharge ICD-9-CM diagnosis codes from insurance claims data as follows: 1) MI or acute coronary syndrome (ICD-9-CM codes 410.xx, 411.1, and 411.8x); 2) ischemic and hemorrhagic stroke (430, 431, 432.x, 433.xx, 434.xx); 3) heart failure (428.xx); or 4) peripheral artery disease (440.21 and 443.9). Here we use the convention 410.xx to denote all codes with category 410 regardless of subcategory and subclassification, and similarly for other categories and subcategories. Because we use insurance claims data, we note that we are able to infer if a patient had a CV event but sought care at an out-of-network hospital. In addition to using diagnosis codes to infer if a CV event occurred, we considered a patient to have experienced a CV event if the cause of death listed on the death certificate included MI or stroke.

5.5. Methods to handle observations in which event status is unknown

Subjects who experienced an event within five years were recorded as E = 1, and those with at least 5 years of event-free follow-up were recorded as E = 0. Subjects who were event-free but censored before accruing 5 years of follow-up have E unknown. We applied and evaluated four variants of each of the machine learning techniques described in Section 3 to our data. The variants differ in their handling of subjects with E unknown:

1. Set E = 0 if E is unknown. Techniques using this strategy are denoted with the suffix -Zero.

2. Discard observations with E unknown. Techniques using this strategy are given the suffix - Discard.

3. Use IPCW on observations with E known. The resulting techniques, as described in Section 3, have the suffix -IPCW.

4. “Split” observations with E unknown into two observations with E = 1 and E = 0 with weights based on marginal survival probability. The resulting techniques, as described subsequently, have the suffix -Split.

The final technique of splitting observations with E unknown was described by Štajduhar and Dalbelo-Bašić [27]. For each observation i in the training set for which E_i is unknown, we create two observations, one with E = 1 and the other with E = 0, but with the same features X_i. Let $\widehat{F}$(t) be the Kaplan-Meier estimator of the survival probability at time t,

where d_j is the number of subjects who are observed to experience the event at time V_j, and n_j is the number of subjects “at risk” for the event (i.e., not yet censored or experienced an event) at time V_j. Then, if E is unknown for instance i in the training set, the weight for the imputed observation with E = 0 is $\widehat{F}$(τ)/$\widehat{F}$(V_i) (an estimate of the conditional probability that E = 0), and the weight for the imputed observation with E = 1 is 1 − $\widehat{F}$(τ)/$\widehat{F}$(V_i). The weights are implemented in the analysis in the same way as the IPC weights. These weights are advantageous because all observations receive non-zero weights and are used in the analysis unlike IPC weights.

5.6. Machine learning methods and implementation details

We now provide some implementation details for the various machine learning techniques. All models were trained on 75% of the sample observations. Code to implement these machine learning methods using each of the four techniques described above to handle observations where the event status is unknown is available as a Github repository from the first author (@docvock).

5.7. Risk prediction evaluation metrics

All measures of predictive performance were assessed on the hold-out test set which was 25% of the sample. To calculate the calibration statistic, we defined five risk strata based on clinically relevant cutoffs for the risk of experiencing a cardiovascular event within 5 years: 0-5%, 5-10%, 10-15%, 15-20% and > 20% [4]. Therefore, a statistical test for the null hypothesis that a model is well-calibrated would reject the null at a 5% significance level if the statistic exceeds 7.81. To calculate the cNRI, we used the same five risk strata as for the calibration statistic. The C-statistic was calculated as described in Section 4.2.

7.1. Limitations

The statistical validity of IPCW rests on several assumptions, in particular that the censoring time is independent of both the event time and patient features. This is a plausible assumption for EHD, where censoring typically occurs for reasons unrelated to a person’s health status, but the assumption is much less plausible in other contexts. For example, if data were collected from a small regional hospital, patients with severe health problems might be censored because they went to a larger facility to seek care. Additionally, IPCW can be inefficient as those subjects for whom the event status is not known are given a weight of zero and do not contribute directly to the estimation of the risk (although those subjects do contribute information to estimate the weights).

Additionally, the example we considered here used a relatively modest number of features to predict CV events as the number of covariates which are currently and consistently measured across all adults in the primary care clinics is relatively small. However, this may change as the complexity, completeness, and quality of EMR data increases.