Regression Trees for Cumulative Incidence Functions

Let $T^{(m)}$ be the time to failure for failure type $m=1,\ldots,K,K\geq 2$ and let $W$ be a vector of $p$ covariates, where $W\in{\cal S}\subset\mathbb{R}^{p}$. Let $T=\min(T^{(1)},\ldots,T^{(K)})$ be the minimum of all latent failure times; it is assumed $(T^{(1)},\ldots,T^{(K)})$ has a bounded multivariate density function and hence that $T$ is continuous. Then, in the absence of other loss to follow-up, $F=(T,W,M)$ is assumed to be the fully observed (or full) data for a subject, where $M$ is the event type corresponding to $T.$ This assumption implies that $(T^{(M)},M,W)$ is observed and, in addition, that $T^{(m)}>T^{(M)}$ for every $m\neq M;$ however, $T^{(m)}$ itself is not assumed to be observed for $m\neq M$. Let $\mathcal{F}=(F_{1},\ldots,F_{n})$ be the full data observed on $n$ independent subjects, where $F_{i}=(T_{i},W_{i},M_{i}),i=1,\ldots,n$ are assumed to be identically distributed (i.i.d.). The dependence structure of $(T^{(1)},\ldots,T^{(K)})$ is left unspecifed.

Let $\Psi_{0}=\{\psi_{0m}(t;w)=P(T\leq t,M=m|W=w),t>0,m=1,\ldots,K\}.$ The set of CIFs $\Psi_{0}$ can be estimated from the data ${\cal F}$ using any of several suitable parametric or semiparametric methods without further assumptions (e.g., independence of $T^{(1)},\ldots,T^{(K)}$). This section introduces a loss-based method for estimating $\psi_{0m}(t;w)$ for a fixed cause $m$ and time point $t>0$ in the case where $\psi_{0m}(t;w)$ is piecewise constant as a function of $W.$ This is a key step in our proposed approach to building a regression tree to estimate $\psi_{0m}(t;w).$

Let $\mathcal{N}_{1},\ldots,\mathcal{N}_{L}$ form a partition of ${\cal S}.$ Define $\beta_{0lm}(t)=P(T\leq t,M=m|W\in\mathcal{N}_{l})$ and suppose $\psi_{0m}(t;w)=\sum_{l=1}^{L}\beta_{0lm}(t)I\{W\in\mathcal{N}_{l}\};$ that is, subjects falling into partition $\mathcal{N}_{l}$ all share the same CIF $\beta_{0lm}(t).$ Define $Z_{m}(t)=I(T\leq t,M=m)$ and let $\psi_{m}(t;w)=\sum_{l=1}^{L}\beta_{lm}(t)I\{W\in\mathcal{N}_{l}\}$ be a model for $\psi_{0m}(t;w).$ Then, fixing both $t>0$ and $m,$ the Brier loss (cf., brier1950verification) $L_{m,t}^{full}(F,\psi_{m})=\{Z_{m}(t)-\psi_{m}(t;w)\}^{2}=\sum_{l=1}^{L}I\{W\in\mathcal{N}_{l}\}\{Z_{m}(t)-\beta_{lm}(t)\}^{2}$ is an unbiased estimator of the risk $\Re(t,\psi_{m})=E[\sum_{l=1}^{L}I\{W\in\mathcal{N}_{l}\}\{Z_{m}(t)-\beta_{lm}(t)\}^{2}],$ or equivalently, $\Re(t,\psi_{m})=\sum_{l=1}^{L}P\{W\in\mathcal{N}_{l}\}\{\beta_{0lm}(t)-\beta_{lm}(t)\}^{2}.$ Assuming that ${\cal F}$ is observed,

is also an unbiased estimator of $\Re(t,\psi_{m}).$ When considered as a function of the set of scalar parameters $\beta_{lm}(t),l=1,\ldots,L$ (i.e., for fixed $m$ and $t$), the empirical Brier loss (1) is minimized when $\psi_{m}(t;w)=\hat{\psi}_{m}(t;w)=\sum_{l=1}^{L}I\{W_{i}\in\mathcal{N}_{l}\}\hat{\beta}_{lm}(t),$ where

is a nonparametric estimate for $\beta_{0lm}(t),m=1,\ldots,K.$

The CART algorithm of breiman1984classification estimates a regression tree as follows:

1. 1.

   Using recursive binary partitioning, grow a maximal tree by selecting a (covariate, cutpoint) combination at every stage that minimizes an appropriate loss function;

2. 2.

   Using cross-validation, select the best tree from the sequence of candidate trees generated by Step 1 via cost complexity pruning (i.e., using penalized loss).

In its most commonly used form, CART estimates the conditional mean response as a piecewise constant function on the covariate space, making all decisions on the basis of minimizing squared error loss. The resulting tree-structured regression function estimates the predicted response in each terminal node using the sample mean of the observations falling into that node. This process generalizes in a straightforward way to more general loss functions.

In the absence of censoring and under a piecewise constant model

for $\psi_{0m}(t;w),$ Section 2.2 shows that a nonparametric estimate for $\psi_{0m}(t;w)$ for a given cause $m$ at a fixed $t>0$ can be obtained by minimizing the loss (1), a problem equivalent to estimating the conditional mean response from the modified dataset ${\cal F}_{red,t}=\{(Z_{im}(t),W^{\prime}_{i})^{\prime},i=1,\ldots n\}$ by minimizing the squared error loss (1). Thus, any implementation of CART for squared error loss (e.g., rpart) applied to ${\cal F}_{red,t}$ will produce a regression tree estimate of $\psi_{0m}(t;w).$ In particular, CART estimates $L$ and the associated terminal nodes $\{\mathcal{N}_{1},\ldots,\mathcal{N}_{L}\}$ from the data ${\cal F}_{red,t},$ and within each terminal node, estimates $\psi_{0m}(t|w)$ by (2). While statistically inefficient, this process can be repeated for each $m=1,\ldots,K$ and any $t>0$ to generate nonparametric estimates of the CIF at a given set of time points for every cause.

Let $C$ be a continuous right-censoring random variable that, given $W$, is statistically independent of $(T,M).$ For a given subject, assume that instead of $F$ we only observe $O=\{\tilde{T},\Delta,M\Delta,W\},$ where $\tilde{T}=\min(T,C)$ and $\Delta=I(T\leq C)$ is the (any) event indicator. The observed data on $n$ i.i.d. subjects is $\mathcal{O}=(O_{1},\ldots,O_{n}).$ Similarly to the case where $K=1,$ random censoring on $T$ permits estimation of the CIF from the data ${\cal O}$.

Fix $t>0$ and let $t^{*}\geq t.$ Define $G_{0}(s|W)=P(C\geq s|W)$ for any $s\geq 0$ and suppose that $G_{0}(T_{i}|W_{i})\geq\epsilon$ almost surely for some $\epsilon>0$. Let $\Delta_{i}(t^{*})=I(T_{i}(t^{*})\leq C)$ and $T_{i}(t^{*})=T_{i}\wedge t^{*};$ in addition, define $\tilde{Z}_{im}(t)=I(\tilde{T}_{i}\leq t,M_{i}=m),$ $i=1,\ldots,n.$ Easy calculations then show

for a fixed $\psi_{m}(t;w).$ This risk equivalence motivates the construction of an IPCW-type observed data loss function that is an unbiased estimator for $\Re(t,\psi_{m}).$ Define for any suitable survivor function $G(\cdot|\cdot)$

then, $L_{m,t}^{ipcw}({\cal O},\psi_{m};t^{*},G_{0})$ has the same risk as (1) and it follows that (4) is minimized by

Observe that (4) and (5) respectively reduce to (1) and (2) when censoring is absent.

When $K=1$ and $t^{*}=\infty,$ we have $\Delta_{i}(\infty)=\Delta_{i}$ and $T_{i}(\infty)=T_{i};$ the loss function (4) is then just a special case of that considered in molinaro2004tree (see also steingrimsson2016doubly). Similarly, for $K=1$ and setting $t^{*}=t,$ the loss function (4) is just that considered in lostritto2012partitioning. Hence, (4) extends several existing loss functions to the problem of estimating a CIF. In practice, an estimator $\hat{G}(\cdot|\cdot)$ for $G_{0}(\cdot|\cdot)$ is used in (4); popular approaches here include product-limit estimators derived from the Kaplan-Meier and Cox regression estimation procedures.

Given $t>0,$ different choices of $t^{*}\geq$ in (4) generate different losses, hence different CART algorithms. However, there are just two important choices of $t^{*}\geq t$ in (4): $t^{*}=\infty$ (standard IPCW; IPCW${}_{\!1}$) and $t^{*}=t$ (modified IPCW; IPCW${}_{\!2}$). In selecting $t^{*}=t,$ any observations that are censored after time $t$ can still contribute (all-cause failure) information to (4); as $t^{*}\rightarrow\infty$, the influence of these observations eventually vanishes. Consequently, the IPCW${}_{\!2}$ loss uses more of the observed data in calculating the loss function and associated minimizers compared to the IPCW${}_{\!1}$ loss and one can therefore expect a regression tree built using (4) for $t^{*}=t$ to perform as well or better than one derived from (4) with $t^{*}=\infty$. The use of $t^{*}=t$ in place of $t^{*}=\infty$ has an additional practical advantage: since $T_{i}(t^{*})\leq T_{i}$ for every $i,$ we have $\hat{G}(\tilde{T}_{i}(t^{*})|W_{i})\geq\hat{G}(\tilde{T}_{i}|W_{i}),$ making it easier to satisfy the empirical positivity condition required for implementation.

Covariates $W_{1},W_{2},W_{3},W_{4},W_{5},W_{6},W_{7},W_{8},W_{9},W_{10}\sim Unif(0,1)$ are generated independently of each other; let $\boldsymbol{W}=(W_{1},W_{2},W_{3},W_{4},W_{5},W_{6},W_{7},W_{8},W_{9},W_{10})^{\prime}.$ The underlying competing risks model assumes $K=2$ and takes the form

where $Z(\boldsymbol{W})=I(W_{1}\leq 0.5\ \&\ W_{2}>0.5)$ and $\beta_{1}$ and $\beta_{2}$ are regression coefficients (cf., fine1999proportional). Note that this formulation satisfies the additivity constraint $\psi_{01}(\infty;\boldsymbol{W})+\psi_{02}(\infty;\boldsymbol{W})=1.$ Three settings are considered: (i) $\beta_{1}=3$ and $\beta_{2}=-0.5$ with $p=0.3$ (High Signal); (ii) $\beta_{1}=2$ and $\beta_{2}=-0.5$ with $p=0.3$ (Medium Signal); and, (iii) $\beta_{10}=1.5$ and $\beta_{20}=-0.5$ with $p=0.3$ (Low Signal). The corresponding cumulative incidence functions for these settings are shown in Figure 1.

Similarly to steingrimsson2016doubly, we respectively generate training sets of size 250, 500 and 1000 (with noninformative censoring) for the purposes of estimation and an independent uncensored test set of size 2000 for evaluating performance. Censoring is exponentially distributed with a rate $\gamma$ chosen to give approximately 50% censoring on $T$. For each training set size, 500 simulations are run for each of the three settings.

For assessing estimator performance, it is assumed there is an underlying true tree structure to be recovered. Let $\hat{\psi}_{1}(t|W_{i})$ be the tree-predicted cumulative incidence for the cause of interest for a subject having covariates $W_{i}$; let $\psi_{01}(t|W_{i})$ be the corresponding true cumulative incidence function. Then, as a measure of predictive performance, we consider

where $n_{test}$ is the number of observations in the test dataset. We also use several other performance measures focused more on the identification of the underlying tree structure:

• •

  $|$fitted size$-$true size$|$ : This quantity measures difference between the size of fitted tree and the size of the true tree, where size denotes the number of terminal nodes. It can be seen from Figure 1 that the true tree size is 3.

• •

  NSP: This quantity is defined as $(n_{sim})^{-1}\sum_{i=1}^{n_{sim}}NS_{i}$ where $NS_{i}$ represents number of times any of the covariates $W_{3},\ldots,W_{10}$ (i.e., covariates that do not affect the true CIF) appear in simulation run $i$.

• •

  PCSP : This quantity is defined as $(n_{sim})^{-1}\sum_{i=1}^{n_{sim}}CT_{i}$ where $CT_{i}=1$ if the fitted tree in simulation $i$ has the same covariates and the same number of splits as the true tree that defines the CIF; otherwise $CT_{i}=0$. While it can happen that a fitted tree may involve only the covariates $W_{1}$ and $W_{2}$, it cannot be expected to have exactly the same splits as the true tree (especially for a continuous covariate); thus, in this measure, a ‘correct tree’ is one that looks like the true tree but may involve different split points.

The IPCW and doubly robust Brier loss functions require estimation of $G_{0}(\cdot|\cdot).$ The Buckley-James and doubly robust Brier loss functions both rely on the function $y_{m}(\cdot;t,w,\Psi)$ in (8), the optimal choice requiring specification of both $\psi_{0m}(u|w)$ and the event-free survival function $\bar{\psi}_{0}(u|w)=P_{\Psi_{0}}(T>u|w)$ for $u>0,w\in{\cal S}$.

Estimation of $G_{0}(\cdot|\cdot)$ is comparatively straightforward; for example, one may use regression procedures (e.g., Cox regression models, random survival forests) or product-limit estimators as appropriate. For all simulations considered here, the censoring distribution $G_{0}(\cdot|\cdot)$ is estimated by the Kaplan-Meier method. For building an IPCW-type tree with $t^{*}=\infty$ in (4), the resulting censoring weight may not satisfy the required positivity assumption. Hence, similarly to steingrimsson2016doubly, we replace $t^{*}=\infty$ by $t^{*}=\tau$ when calculating IPCW${}_{\!1},$ where $\hat{G}(\tilde{T}_{i}(\tau))\geq 0.05,i=1,\ldots,n$ (i.e., the marginal truncation rate of observed survival times is 5%). No such modification is required when calculating IPCW${}_{\!2},$ that is, when computing (4) with $t^{*}=t=t_{j},j=1,2,3.$

Calculation of the Buckley-James and doubly robust Brier loss functions in practice requires using an estimated model $\hat{\Psi}$ in place of $\Psi_{0}$ when computing (8). Parametric models have been proposed that ensure the compatibility between the models used for $\psi_{m}(\cdot|\cdot)$ and $\bar{\psi}(\cdot|\cdot)$ as well as the ability to calculate probabilities for every $u>0;$ see, for example, jeong2006direct, jeong2007parametric, cheng2009modeling, and shi2013constrained. Compared to existing semiparametric methods (e.g., fine1999proportional; scheike2008predicting; eriksson2015proportional), the parametric likelihood methods of jeong2006direct and jeong2007parametric enable researchers to estimate two or more CIFs jointly. However, cheng2009modeling and shi2013constrained show that even these methods fail to enforce the additivity constraint $\sum_{m=1}^{M}P_{\Psi}(T\leq\infty,M=m;W=w)=1$ and propose methods that resolve this inconsistency.

Since our focus is on a particular cause $m,$ we will without loss of generality assume that $m=1$ and that $K=2$ (i.e., $m=2$ captures all causes $m\neq 1$); in this case, specification of $\psi_{1}(\cdot|\cdot)$ and $\bar{\psi}(\cdot|\cdot)$ is equivalent to specifying $\Psi=(\psi_{1}(\cdot|\cdot),\psi_{2}(\cdot|\cdot)).$ Our simulation study considers the following estimators $\hat{\Psi}$ when $K=2:$

1. 1.

   RF(npar) : random survival forests as proposed in ishwaran2014random.

2. 2.

   RF+lgtc: a hybrid method of random survival forests and parametric modeling;

3. 3.

   FG(true): the simulation model of Section 4.1 is an example of the class of models considered in jeong2007parametric. For $m=1,2,$ $\psi_{m}(\cdot|\cdot)$ are estimated by maximum likelihood under a correctly specified parametric model.

4. 4.

   lgtc+godds: for $m=1,2,$ $\psi_{m}(\cdot|\cdot)$ is estimated by maximum likelihood using the parametric cumulative incidence model proposed in shi2013constrained.

Further details on the RF(npar) and RF+lgtc estimators, including the parametric models used in #2 and #4 above, may be found in Section A.3 of the Appendix. We anticipate that FG(true) will lead to the best performance because the underlying models for the event of interest and for censoring are both specified correctly. As in steingrimsson2016doubly, the estimates $\hat{G}(\cdot|\cdot)$ and $\hat{\Psi}$ depend only on ${\cal O}$ and are computed in advance of the process of building the desired regression tree for estimating $\psi_{01}(\cdot|\cdot).$

Each combination of loss function and method for estimating $G_{0}$ and/or $\Psi_{0}$ leads to a distinct algorithm. In this section, we show results for the IPCW-type estimators with $t^{*}=t_{j},j=1,2,3$ and $t^{*}=\tau$; the Buckley James Brier loss, where RF(npar) (algorithm BJ-RF (npar)) and FG(true) (algorithm BJ-FG (true)) are respectively used to estimate $\Psi_{0}$; and, the doubly robust Brier loss function, where RF(npar) (algorithm DR-RF (npar)) and FG(true) (algorithm DR-FG (true)) are respectively used to estimate $\Psi_{0}$. Results using RF+lgtc and lgtc+godds for estimating $\Psi_{0}$ are summarized in the Appendix (Section A.5).

In fitting the trees, we consider three fixed time points $t_{1}$, $t_{2}$ and $t_{3}$, respectively representing the 25th, 50th and 75th quantile of the true marginal distribution of $T$. These times are used for building individual trees as well as for composite loss functions with the respective training set sizes $n=250,500$ and $1000$. Below, we focus on the results with $n=500$ using composite loss functions as described in Section 3.3.2. Table 1 summarizes the tree fitting performance measures for estimating $\psi_{01}(\cdot|\cdot)$ in (13) for $n=500$ using a composite loss function having weights $w_{1}=w_{2}=w_{3}=1/3$. For the chosen metrics, and irrespective of signal strength, there is a clear benefit to using augmented loss functions over IPCW loss functions, the exception being BJ-RF (npar) (i.e., the Buckley-James loss with $\hat{\Psi}$ estimated by random survival forests). Overall, the BJ-FG (true) algorithm that uses the correct model for $\Psi_{0}$ performs best, followed by DR-FG (true), DR-RF (npar), IPCW${}_{\!2},$ BJ-RF (npar) and then IPCW${}_{\!1}$. Prediction error performance is considered in Figure 2, which shows boxplots of the mean squared prediction error given in (15) under equally weighted composite loss. Here, the BJ-FG (true) algorithm again performs uniformly best, with the augmented methods DR-FG (true), DR-RF (npar), and BJ-RF (npar) all performing similarly, followed by IPCW${}_{\!2}$ and then IPCW${}_{\!1}.$

In general, the methods that incorporate information beyond that used by IPCW₁ exhibit significant improvement in performance, with the methodology for estimating quantities needed to compute the augmented loss function playing a smaller role, particularly as sample size increases. Considering only methods where no knowledge of $\Psi_{0}$ is used, the doubly robust loss function also provides gains over both the Buckley-James and IPCW${}_{\!2}$ loss functions; however, the gains achieved are substantially less compared to the gains over IPCW${}_{\!1}$. This phenomenon is expected and can be explained by noting that (i) the censoring distribution is being consistently estimated in all cases; (ii) IPCW${}_{\!2}$ can be viewed as an augmented version of IPCW${}_{\!1};$ and, (iii) IPCW${}_{\!2}$ only requires estimating one infinite dimensional parameter.