Dynamic Risk Prediction Triggered by Intermediate Events Using Survival Tree Ensembles

Cystic fibrosis (CF) is a genetic disease characterized by a progressive, irreversible decline in lung function caused by chronic microbial infections of the airways. Despite recent advances in diagnosis and treatment, the burden of CF care remains high, and most patients succumb to respiratory failure. There is currently no cure for CF, so early prevention of lung disease for high-risk patients are essential for successful disease management. The goal of this research is to develop flexible and accurate event risk prediction algorithms for abnormal lung function in pediatric CF patients by exploiting the rich longitudinal information made available by the Cystic Fibrosis Foundation Patient Registry (CFFPR).

The CFFPR is a large electronic health record database that collects encounter-based records of over 300 unique variables on patients from over 120 accredited CF care centers in the United States (Knapp et al., 2016). The CFFPR contains detailed information on potential risk factors, including symptoms, pulmonary infections, medications, test results, and medical history. Analyses of CFFPR suggested that the variability in spirometry measurements over time is highly predictive of subsequent lung function decline (Morgan et al., 2016). Moreover, the acquisition of chronic, mucoid, or multidrug-resistant subtypes of Pseudomonas aeruginosa (PA) leads to more severe pulmonary disease, accelerating the decline in lung function (McGarry et al., 2020). In this paper, the event of interest is the progressive loss of lung function, defined as the first time that the percent predicted forced expiratory volume in 1 second (ppFEV1) drops below 80% in CFFPR. Since risk factors such as weight and height in pediatric patients can change substantially over time, models with baseline predictors have limited potential for long-term prognosis. Incorporating repeated measurements and intermediate clinical events would reflect ongoing CF management and result in more accurate prediction.

To incorporate the longitudinal patient information in risk prediction, one major approach is joint modeling (see, for example, Rizopoulos, 2011; Taylor et al., 2013). Under the joint modeling framework, a longitudinal submodel for the time-dependent variables and a survival submodel for the time-to-event outcome are postulated; the sub-models are typically linked via latent variables. Such a model formulation provides a complete specification of the joint distribution, based on which the survival probability given the history of longitudinal measurements can be derived. Most joint modeling methods consider a single, continuous time-dependent variable. Although attempts have been made to incorporate multiple time-dependent predictor variables (Proust-Lima et al., 2016; Wang et al., 2017), correct specification of the model forms for all the time-dependent covariates and their associations with the event outcome remains a major challenge. Moreover, it is not clear how existing joint modeling approaches can further incorporate the information on the multiple intermediate events, such as the acquisition of different subtypes of PA, in risk prediction.

Another major approach that can account for longitudinal predictors is landmark analysis, where models are constructed at pre-specified landmark times to predict the event risk in a future time interval. For example, at each landmark time, one may postulate a working Cox model with appropriate summaries of the covariate history up to the landmark time (e.g., last observed values) as predictors and then fit the Cox model using data from subjects who are at risk of the event. The estimation can either be performed using a separate model at each landmark time point or a supermodel for all landmark time points (van Houwelingen, 2007; van Houwelingen and Putter, 2008, 2011). This way, multiple and mixed type time-dependent predictors can be easily incorporated. Moreover, to better exploit the repeated measurements and to handle measurement errors, one may also consider a two-stage landmark approach (Rizopoulos et al., 2017; Sweeting et al., 2017; Ferrer et al., 2019): in the first step, mixed-effects models are used to model the longitudinal predictors; in the second step, functions of the best linear unbiased prediction (BLUP) estimator of the random effects are included as predictors of the landmark Cox model. Other than Cox models, Parast et al. (2012) considered time-varying coefficient models to incorporate a single intermediate event and multiple biomarker measurements. Zheng and Heagerty (2005), Maziarz et al. (2017), and Zhu et al. (2019) further considered the impact of informative observation times of repeated measurements on future risk.

Direct application of the existing landmark analysis method to the CFFPR data may not be ideal for the following reasons: first, imposing semiparametric working models at different landmark times may result in incompatible models and inconsistent predictions (Jewell and Nielsen, 1993; Rizopoulos et al., 2017). In other words, a joint distribution of predictors and event times that satisfies the models at all the landmark times simultaneously may not exist. Second, the specification of how the predictor history affects the future event risk may require deep clinical insight. For example, researchers have shown that various summaries of the repeated measurements, including the variability (Morgan et al., 2016), the rate of change (Mannino et al., 2006), and the area under the trajectory curve (Domanski et al., 2020), can serve as important predictors of disease risks, while the last observed value has been commonly used in the statistical literature. Therefore, nonparametric statistical learning methods are appealing in landmark prediction, because they require minimal model assumptions and have the potential to deal with a large number of complicated predictors. Tanner et al. (2021) applied super learners for landmark prediction in CF patients, where discrete time survival analysis were conducted via the use of machine learning algorithms for binary outcomes.

In this paper, we propose a unified framework for landmark prediction using survival tree ensembles, where the landmark times can be subject-specific. A subject-specific landmark can be defined by an intermediate clinical event that modifies patients’ risk profiles and triggers the need for an updated evaluation of future risk. In our application, the acquisition of chronic PA usually leads to accelerated deterioration in the pulmonary function and serves as a natural landmark. When the landmark time is random, the number of observed predictors at the landmark time often varies across subjects, creating analytical challenges in fully utilizing the available information. Moreover, unlike static risk prediction models where baseline predictors are completely observed, the observation of the time-dependent predictors is further subject to right censoring. To tackle these problems, we propose a risk-set-based approach to handle the possibly censored predictors. To avoid the instability issue of a single tree, we propose a novel ensemble procedure based on averaging unbiased martingale estimating equations derived from individual trees. Our ensemble method is different from existing ensemble methods that directly average the cumulative hazard predictions and has strong empirical performances in dealing with censored data.

The rest of this article is organized as follows. In Section 2, we introduce a landmark prediction framework that incorporates the repeated measurements and intermediate events. In Section 3, we propose tree-based ensemble methods to deal with censored predictors and outcomes. We propose a concordance measure to evaluate the prediction performance in Section 4 and define a permutation variable importance measure in Section 5. The proposed methods are evaluated by extensive simulation studies in Section 6 and are applied to the CFFPR data in Section 7. We conclude the paper with a discussion in Section 8.

In contrast to static risk prediction methods that output a conditional survival function given baseline predictors, dynamic landmark prediction focuses on the survival function conditioning on the predictor history up to the landmark time. Since history information involves complicated stochastic processes, challenges arise as to how to partition the history processes when applying tree-based methods. In what follows, we first introduce a generalized definition of the landmark survival function, starting from either a fixed or subject-specific landmark time. We then express the history information as a fixed-length predictor vector on which recursive partition can be applied.

Denote by $T$ a continuous failure event time and by ${T_{\rm L}}$ a landmark time. The landmark is selected a priori and is usually clinically meaningful. We allow ${T_{\rm L}}$ to be either fixed or subject-specific. We focus on the subpopulation that is free of the failure event at ${T_{\rm L}}$ and predict the risk after ${T_{\rm L}}$. Denote by ${\bm{Z}}$ the baseline predictors and denote by ${\mathcal{H}}(t)$ other information observed on $[0,t]$. Our goal is to predict the probability conditioning on all the available information up to ${T_{\rm L}}$, that is,

To illustrate the observed history ${\mathcal{H}}(t)$, we consider two types of predictors that are available in the CFFPR data. The first type of predictors is repeated measurements of time-dependent variables such as weight and ppFEV1. It is worthwhile to point out that both internal and external time-dependent predictors can be included, as only their history up to ${T_{\rm L}}$ will be used. We denote this type of predictors by ${\bm{W}}(t)$, a $q$-dimensional vector of time-dependent variables, and assume ${\bm{W}}(\cdot)$ is available at fixed time points $t_{1},t_{2},\ldots,t_{K}$. The observed history up to $t$ is ${\mathcal{H}}_{W}(t)=\{{\bm{W}}(s)dO(s),0<s\leq t\}$, where $O(t)$ is a counting process that jumps by one when ${\bm{W}}(\cdot)$ is measured (i.e., $dO(t_{k})=1$ for $k=1,\ldots,K$). The second type of predictors is the timings of intermediate clinical events such as pseudomonas infections. Denote by $U_{j}$ the time to the $j$th intermediate event, $j=1,\ldots,J$. The observed history up to $t$ is ${\mathcal{H}}_{U}(t)=\{I(U_{j}\leq s),0<s\leq t,j=1,\ldots,J\}$. Collectively, we have a system of history processes ${\mathcal{H}}(t)=({\mathcal{H}}_{W}(t),{\mathcal{H}}_{U}(t))$.

In our framework, both $t_{k}$ and $U_{j}$ can serve as landmark times. Due to the stochastic nature of $U_{j}$, the order of $\{t_{k},U_{j},k=1,\ldots,K,j=1,\ldots,J\}$ can not be pre-determined. As a result, the number of available predictors at a given landmark time can vary across subjects. For illustration, we consider one fixed time point $t_{1}$ at age 7 and one intermediate event chronic PA (cPA), of which the occurrence time is denoted by $U_{1}$. Figure 1 depicts the observed data of two study subjects. At $t_{1}=7$, subject 1 has experienced cPA ($U_{1}=4.8\leq t_{1}$), while subject 2 remains free of cPA ($U_{1}=10.2>t_{1}$). Let ${\bm{W}}(t)$ be the body weight measured at time $t$ (i.e., $q=1$). The probabilities of interest are given as follows:

1. (I)

   At a fixed landmark time ${T_{\rm L}}=t_{1}$, we predict the risk at $t_{1}+t$, $t>0$, among those who are at risk, that is, $T\geq t_{1}$. Note that subjects in the risk set may or may not have experienced the intermediate event prior to $t_{1}$. Given ${\bm{W}}(t_{1})$ and the partially observed $U_{1}$, the conditional survival probability (1) can be reexpressed as

   $\displaystyle\begin{cases}P(T\geq t+t_{1}\mid T\geq t_{1},{\bm{Z}},{\bm{W}}(t_{1}),U_{1}),&\text{if~{}}U_{1}\leq t_{1},\\ P(T\geq t+t_{1}\mid T\geq t_{1},{\bm{Z}},{\bm{W}}(t_{1}),U_{1}>t_{1}),&\text{otherwise}.\end{cases}$

   In other words, at $t_{1}$, we output the former for subjects who experience the intermediate event prior to $t_{1}$ (subject 1, Figure 1(a)), while output the latter for others (subject 2, Figure 1(b)).

2. (II)

   At a random landmark time ${T_{\rm L}}=U_{1}$, we predict the risk for subjects who have experienced the intermediate event and are free of the failure event. The predictor value ${\bm{W}}(t_{1})$ is available only if $U_{1}\geq t_{1}$. In this case, we predict

   $\displaystyle\begin{cases}P(T\geq t+U_{1}\mid T\geq U_{1},{\bm{Z}},U_{1}),&\text{if~{}}U_{1}\leq t_{1},\\ P(T\geq t+U_{1}\mid T\geq U_{1},{\bm{Z}},{\bm{W}}(t_{1}),U_{1}),&\text{otherwise}.\end{cases}$

   Therefore, at $U_{1}$, we output the former for subjects whose ${\bm{W}}(t_{1})$ is observed after $U_{1}$ (Subject 1, Figure 1(c)), while output the latter for others (Subject 2, Figure 1(d)).

In the above example, the observed predictors vary across subjects. We then represent the history ${\mathcal{H}}({T_{\rm L}})$ as a vector with a fixed length, so that tree-based methods can be applied to estimate the probability in (1). Define the complete predictors ${\bm{X}}=\{{\bm{W}}(t_{1}),\ldots,{\bm{W}}(t_{K}),U_{1},\ldots,U_{J}\}$. The information in ${\bm{X}}$ may not be fully available at a given landmark time. We define the available information up to $t$ by ${\bm{X}}(t)=\{{\bm{W}}(t_{1},t),\ldots,$ ${\bm{W}}(t_{K},t),U_{1}(t),\ldots,U_{J}(t)\}$, where

By a slight abuse of notation, we write $U_{j}(t)=t^{+}$ if $U_{j}>t$ and ${\bm{W}}(t_{k},t)=\textbf{NA}_{q}$ if $t_{k}>t$, with $\textbf{NA}_{q}\overset{\rm def}{=}(\rm NA,\ldots,NA)$ denoting a non-numeric $q$-dimensional vector. Here an NA value indicates that the covariate value is collected after the landmark time and thus should not be used for constructing prediction models. In other words, under our setting NA is treated as an attribute rather than missing data, as the target probability is not conditional on ${\bm{W}}(t_{k})$ for $t_{k}>{T_{\rm L}}$. The covariate history ${\mathcal{H}}({T_{\rm L}})$ can then be expressed as ${\bm{X}}({T_{\rm L}})$, which is a $(qK+J)$-dimensional vector. This way, a covariate not being observed is predictive of the outcome, and the target survival probability function can be expressed as follows:

In the example depicted in Figure 1, we have ${\bm{X}}(t)=\{{\bm{W}}(t_{1},t),U_{1}(t)\}$. For ${T_{\rm L}}=t_{1}$, the predictor values in Figure 1(a) and 1(b) correspond to ${\bm{x}}=(19.7,4.8)$ and ${\bm{x}}=(21.5,7^{+})$, respectively. For ${T_{\rm L}}=U_{1}$, the predictor values in Figures 1(c) and 1(d) correspond to ${\bm{x}}=({\rm NA},4.8)$ and ${\bm{x}}=(21.5,10.2)$, respectively. Since the information at ${T_{\rm L}}$ involves left-bounded intervals and non-numeric values, applying semiparametric methods to estimate $S(t\mid a,{\bm{z}},{\bm{x}})$ is challenging. Hence we propose tree-based methods to handle partially observed predictors.

At a given landmark time, we build a tree-based model to predict future event risk. Survival trees are popular nonparametric tools for risk prediction. The original survival trees take baseline predictors as input variables and output the survival probability conditioning on the baseline predictors (see, for example, Gordon and Olshen, 1985; Ciampi et al., 1986; Segal, 1988; Davis and Anderson, 1989; LeBlanc and Crowley, 1992, 1993; Zhang, 1995; Molinaro et al., 2004; Steingrimsson et al., 2016), and ensemble methods have been applied to address the instability issue of a single tree (Hothorn et al., 2004, 2006; Ishwaran et al., 2008; Zhu and Kosorok, 2012; Steingrimsson et al., 2019). However, existing methods may not be directly applied, because the predictors in ${\bm{X}}$ are not completely observed at ${T_{\rm L}}$, and the available predictors $({T_{\rm L}},{\bm{X}}({T_{\rm L}}))$ are subject to right censoring. In the absence of censoring, we introduce a partition scheme for subjects who are event-free at the landmark time in Section 3.1. To handle censored data, we propose risk-set methods to estimate the partition-based landmark survival probability in Section 3.2 and propose an ensemble procedure in Section 3.3.

To evaluate the performance of the predicted risk score, we extend the cumulative/dynamic receiver operating characteristics (ROC) curves (Heagerty et al., 2000), which has been commonly used when a risk score is based on baseline predictors. We note that ROC and concordance indices for dynamic prediction at a fixed landmark time has been studied in the literature (Rizopoulos et al., 2017; Wang et al., 2017). Here we consider a more general case where the landmark time ${T_{\rm L}}$ is random and subject to right censoring. For $t>0$, subjects with $0\leq T-{T_{\rm L}}<t$ are considered as cases and subjects with $T-{T_{\rm L}}\geq t$ are considered as controls. The ROC curve then evaluates the performance of a risk score that discriminates between subjects who have experienced the events prior to ${T_{\rm L}}+t$ and those who do not.

Let $g(t,{T_{\rm L}},{\bm{Z}},{\bm{X}}({T_{\rm L}}))$ denote a risk score based on $({T_{\rm L}},{\bm{Z}},{\bm{X}}({T_{\rm L}}))$, with a larger value indicating a higher chance of being a case. For each $t>0$, The true positive rate and false positive rate at a threshold of $c$ are defined as follows,

where $\tau_{0}$ is a pre-specified constant. The ROC curve at $t$ is defined as ${\rm ROC}_{t}(p)={\rm TPR}_{t}({\rm FPR}_{t}^{-1}(p))$. Following the arguments of McIntosh and Pepe (2002), it can be shown that $g(t,a,{\bm{z}},{\bm{x}})=1-S(t\mid a,{\bm{z}},{\bm{x}})$ yields the highest ROC curve, which justifies the use of the proposed time-dependent ROC curve. Moreover, the area under the ROC curve is equivalent to the following concordance measure,

where $({T_{\rm L}}_{1},{\bm{Z}}_{1},{\bm{X}}_{1}({T_{\rm L}}_{1}),T_{1})$ and $({T_{\rm L}}_{2},{\bm{Z}}_{2},{\bm{X}}_{2}({T_{\rm L}}_{2}),T_{2})$ are independent pairs of observations, and the second term accounts for potential ties in the risk score.

In practice, one usually builds the model on a training dataset and evaluates its performance on an independent test dataset that are also subject to right censoring. To simplify notation here, we construct the estimator for ${\rm CON}_{t}(g)$ using the observed data introduced in Section 3.2, although ${\rm CON}_{t}$ evaluated using the test data should be used in real applications. Define $d_{ij}(t)=g(t,{Y_{\rm L}}_{j},{\bm{Z}}_{j},{\bm{X}}_{j}({Y_{\rm L}}_{j})){\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-}g(t,{Y_{\rm L}}_{i},{\bm{Z}}_{i},{\bm{X}}_{i}({Y_{\rm L}}_{i}))$. The ${\rm CON}_{t}(g)$ measure can be estimated by

where $\widehat{S}_{C}(t\mid{\bm{z}})$ is an estimator for the conditional censoring distribution $S_{C}(t\mid{\bm{z}})=P(C\geq t\mid{\bm{Z}}={\bm{z}})$. For example, survival trees or forests can be applied to estimate $S_{C}(t\mid{\bm{z}})$; when censoring is completely random, the Kaplan-Meier estimator can also be applied. Under regularity conditions, we show that $\widehat{{\rm CON}}_{t}(g)$ consistently estimates ${{\rm CON}}_{t}(g)$ in the Supplementary Material. In practice, one can either use the concordance at a given time point $t$ or an integrated measure to evaluate the overall concordance on a given time interval $[t_{\rm L},t_{\rm U}]$. In the latter case, a weighted average of concordance on a grid of time points in $[t_{\rm L},t_{\rm U}]$ can be reported. For example, one may assign equal weights to all the time points. Another example is a weight proportional to the denominator in $\widehat{\rm CON}_{t}(g)$, that is, $\sum_{i\neq j}{\Delta_{j}I(0\leq Y_{j}-{Y_{\rm L}}_{j}\leq t<Y_{i}-{Y_{\rm L}}_{i},{{Y_{\rm L}}_{i}\leq\tau_{0},{Y_{\rm L}}_{j}\leq\tau_{0}})}/{{\widehat{S}}_{C}(Y_{j}\mid{\bm{Z}}_{j}){\widehat{S}}_{C}({Y_{\rm L}}_{i}+t\mid{\bm{Z}}_{i})}$, to avoid potentially unstable estimation for very small or large time points.