Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
94,587 characters · 23 sections · 60 citation commands
Censored Quantile Regression Forests
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 {
\affil{Department of Mathematics \\ University of California San Diego}
} \fi
\if10 {
} \fi
{\it Keywords:} Random Forest; Censored quantile regression; Nonparametric regression; Kaplan-Meier estimation. \spacingset{1.45}
In many applications, we want to predict and estimate the effect of a covariate on survival time of interests. Examples include treatment, surgical procedure, or immunization on survival time of patients, who for example, could be individuals who have metastatic breast cancer, military casualties suffering from various injuries, or survival time of infectious diseases. Classically, most datasets have been too small to meaningfully examine the heterogeneity of the data beyond dividing them into a few sub-populations. In the past few years, however, there has been an explosion of experimental settings where it is potentially feasible to explore heterogeneity to its full extent.
An impediment to exploring heterogeneous effects is the fear that scientists with two opposite agendas could hypothetically string together two opposite but coherent results by searching through many different possible models and then reporting only the very extreme ones -- highlighting solely spurious results olken2015promises. Thus, protocols for clinical trials must specify in advance the pre-analysis plans and then learn from the data. However, such restrictions can make it challenging to discover unexpected effects due to heterogeneity. Here, we aim to address this challenge by developing a robust, non-parametric method for estimation in regression settings with censored response variable, which yields consistent estimator that adapts to the heterogeneity of the data and hence can be broadly applied to many different models and achieves further improvement. One example is the accelerated failure time model.
Classical approaches to accelerated failure time model include non-parametric maximum likelihood method, semi-parametric approaches, as well as traditional parametric approaches; see e.g., zeng2007efficient, robins1992semiparametric, robins1992estimation, as well as koul1981regression, wei1992accelerated, huang2007least etc. These methods perform well in settings where the model error is correctly specified, but quickly break down when the distribution of the error has any heterogeneous, asymmetric, or outlier structure. Other non-parametric methods, like doi:10.1093/biomet/68.2.381, doi:10.1093/biomet/85.4.809, and rank-based methods, like jin2003rank, strive to achieve certain level of assumption-lean modeling of the mean in the accelerated failure time model. In this paper, we explore the use of ideas from the machine learning literature to improve the performance of these classical methods in a non-parametric fashion that is adaptive to heterogeneous or heavy-tailed error distributions.
Random forest algorithms introduced by breiman2001random allow for flexible modeling of covariate interactions, and are related to kernels and nearest-neighbor methods in the sense that they make predictions using a weighted average of “nearby” observations. However, random forests depart from the above principles in that they have a data-driven way to determine which nearby observations receive more weight, something that is especially important in environments with many covariates or complex interactions among covariates.
Despite their wide-spread success at estimation and prediction, application of random forests to censored regression models is far from being easily understood. Not all response variables being observed makes it difficult to understand how to evaluate the prediction arising from simple tree structures. Namely, a simple average is only enough for conditional mean estimation without censoring. Median estimates, and quantiles more generally, are not easy to construct. In particular, quantile random forests meinshausen2006quantile have mainly been undeveloped for censored observations.
This paper addresses these limitations, developing a forest-based method for estimation of quantiles for randomly censored observations (right or left censored). More formally, let $T$ be a real-valued latent variable (e.g., survival time) and $X$ be a (possibly high-dimensional) predictor variable. Let $C$ denote the censoring variable, which prevents us from observing all of the information regarding $T$. In left-censored data, we only observe $Y_i = \max(T_i,C_i)$, whereas in right-censored data, we only observe $Y_i = \min(T_i,C_i)$. Our goal is to estimate the $\tau$-th quantile of $T_i$, non-parametrically, using observations $(Y_i,X_i, \delta_i)$ where $\delta_i =\operatorname*{\mathbbm{1}} \{ T_i \leq C_i\}$.
We take on the perspective of random forests as that of an nonparametric, adaptive kernel smoothing method. This interpretation follows work by athey2016generalized, Bloniarz2016SupervisedNF, hothorn2004bagging, li2017forest, and meinshausen2006quantile, and supplements the customary view of forests as an ensemble method (i.e., an average of predictions made by separate trees). In this view, random forest predictions, of the mean for example, evaluated at a new test point $x$, can be represented by \[ \sum_{i=1}^n w_i(X_i,x) Y_i, \qquad \sum_{i=1}^n w_i(X_i,x) =1, \qquad w_i(X_i,x)\geq 0 \] where the weights encode the similarity between the new test point $x$ and the observed covariates $X_i$. It is worth pointing that for the conditional mean, the averaging and weighting views of random forests are equivalent; however, once we move to more general settings, the weighting-based perspective proves substantially more powerful. The goal is then to utilize these local neighborhood weights and quantile regression adjustments to design a new non-parametric quantile estimate of
based on observations $(Y_i, X_i, \delta_i)_{i=1}^n$. In their simplest form, censored quantile forests just take the forest weights, $w_i(X_i,x)$, and use them for quantile regression: \[ \hat q_\tau (x) = \arg\min_q \left\{ \Bigl\| \sum_{i=1}^n w_i(X_i,x) {\mathcal{S}}_\tau (X_i, Y_i, q,x) \Bigl\|_l\right\}, \] where $\| \cdot \|_l$ denotes any $l$-norm with $l\geq 1$. Here, we define a new score function, ${\mathcal{S}}_\tau $, that is censoring and quantile adaptive: \[ {\mathcal{S}}_\tau (X_i, Y_i, q,x)= (1-\tau)\hat{G}(q|x) - \operatorname*{\mathbbm{1}} (Y_i > q) \] with $\hat G$ denoting an estimate of $G (u) = \mathbb{P}(C \geq u| X =x)$, the survival function of the censoring time $C $ given $X$. In the rest of the document we focus on the case of $l=1$; however, results remain true for general cases of $l$ as well.
Formally, we study the performance of the above non-parametric quantile estimator $\hat q_\tau(x)$ of (ref) while considering right-censored observations. All of our results can be trivially extended to the left-censored case. Most of the theoretical work focuses on establishing Theorem (ref) that ensures the consistency of the censored quantile forest estimator $\hat q_\tau(x)$ at a given test point $x$ while allowing minimal assumptions on the regression model as well as the distribution of the model error. This result also allows us to construct prediction intervals regarding (ref) that adapt to the model structure of the data generating process. Censored quantile forest estimator improves on prediction error compared with many state-of-the-art censored regression methods, and yet it retains the flexibility of random forest methods designed for non-censored observations. Finally, we note that our method can also be seen as an improvement over classical non-parametric approaches for censored observations. While the latter only perform well in low-dimensional problems, ours performs well even with a large number of covariates. One important reason is that random forest weights adapt reasonably well to the dimensionality increase, whereas kernel methods suffer from the curse of dimensionality.
There has been a long-time understanding that proportional hazard models, and cox's model, in particular, are especially powerful for right-censored observations and regression problems. However, they do not adapt to possibly left-censored observations; besides, they heavily rely on the proportionality assumption which can sometimes be inappropriate, necessitating stratification of the baseline hazard or some other weakening of the proportional hazards condition koenker2008censored.
A more flexible approach for random censoring problem is to model the conditional quantiles of the response variable directly. This approach offers greater flexibility as it does not restrict the structure of the hazard function koenker2008censored and merely is more intuitive. To estimate the conditional quantiles, portnoy2003censored proposed a recursive method which estimates a sequence of linear conditional quantile functions recursively. It can be treated as a generalization to regression of the Kaplan Meier estimator. Another closely related quantile regression model proposed by peng2008survival instead makes the linkage to the Nelson-Aalen estimator of the cumulative hazard function, upon which they developed a complete asymptotic theory. The closest work to ours is that of doi:10.1080/01621459.2018.1469996 where authors propose a new censored loss function. However, they only discuss the properties of quantile estimates in the case of linear quantile model.
The parametric methods, including those mentioned above, always rely on the linearity assumption on the conditional quantiles, that is,
Here, the $\log$ transformation is arbitrary but popular in survival analysis and can be replaced by any monotone function. This linearity assumption is too restrictive in many cases, especially when data lie on a complex manifold. Therefore, non-parametric methods are necessary and play an important role in modeling data heterogeneity.
In the case of right censoring, most non-parametric recursive partitioning algorithms in the existing literature rely on survival tree or its ensembles. ishwaran2008random proposed random survival forest (RSF) algorithm in which each tree is built by maximizing the between-node log-rank statistic. However, it is not directly estimating the conditional quantiles but instead estimating the cumulative hazard. zhu2012recursively proposed the recursively imputed survival trees (RIST) algorithm with the same splitting criterion for each individual tree but different ensemble scheme. Other similar methods relying on different kinds of survival trees were proposed in gordon1985tree, segal1988regression, davis1989exponential, leblanc1992relative, and leblanc1993survival. All these methods as mentioned above use splitting rules specifically designed to deal with the right censored data. Despite different splitting strategies, they all rely on the proportional hazard assumption and cannot reduce to a loss-based method that might ordinarily be used in the situation with no censoring.
molinaro2004tree proposed a tree method based on the inverse probability censoring robins1994estimation weighted (IPCW) loss function which reduces to the full data loss function used by CART in the absence of censoring. hothorn2005survival then extended the IPCW idea and proposed a forest-type method in which each tree is trained on resampled observations according to inverse probability censoring weights. However, the censored data always get weights zero and hence only uncensored observations will be resampled. As pointed out by robins1994estimation, the inverse probability weighted estimators are inefficient because of their failure to utilize all the information available on observations with missing or partially missing data.
This work aims to build a non-parametric conditional quantile estimator for randomly censored data that reduces to the ordinal quantile forest estimators when full data are observed, and efficiently utilizes all available information on both uncensored and censored observations. Furthermore, it does not require the specific modification of ordinal regression tree (e.g., CART) to survival tree, and hence works on both left and right censored problems.
Fundamentally different from the aforementioned forests methods, in which the censoring information is considered directly in the tree constructing process, our method avoids this complexity and only requires building ordinal regression trees (e.g., CART) for the first step, treating all observations equally. The censoring effects are then taken care of in the second step by solving a locally weighted estimating equation. These local weights can be directly calculated from the random forest constructed in the first step; weights derived from the fraction of trees in which an observation appears in the same leaf as the target value of the covariate vector. This locally weighted view of random forests was previously advocated by hothorn2005survival and Bloniarz2016SupervisedNF; original random forest algorithm breiman2001random utilized ensemble learning literature. Differently from kernel weights, typically employed in local maximum likelihood method, for example, and whose performance suffers greatly whenever the dimensionality of the covariate space is more than two or three, our random-forests weights adapt well to moderate dimensionality increases.
Additional challenges arise due to the random censoring nature of the observations. For fixed censoring, one observes all the censoring values and hence can straightforwardly modify the objective used in the general framework of athey2016generalized, for example. However, it is unclear how to develop a non-parametric estimator that adapts to unknown censoring in the observations.
The paper is organized as follows. In Section (ref) we provide local adaptive nature of random forests and regression adjustment utilizing those weights. In Sections (ref) and (ref), we showcase the development of a new loss function and its power in predicting any conditional quantile of the latent variable $T$ by solving an ingenious estimating equation, which is designed to correct the censoring effect. The Algorithm is then described in details in Section (ref). In Section (ref), we analyze the time complexity of our algorithm and prove consistency of the proposed estimator. Section (ref) contains extensive numerical studies where we compare our algorithm with other forest algorithms on simulated and real censored data sets. Proofs are collected in the Appendix.
The quantile random forest cannot be directly applied to censored data $\{(X_i, Y_i)\}$ because the conditional quantile of $Y$ is different from that of the latent variable $T$ due to the censoring. Moreover, there is no explicitly defined quantile loss function for randomly censored data. In this section, we design a new approach to achieve both tasks.
Let $\theta$ denote the random parameter determining how a tree is grown, and $\{(X_i, Y_i): i=1,\ldots,n\} \in \mathcal{X} \times \mathcal{Y} \subset \mathbb{R}^p \times \mathbb{R}$ denote the training data. For each tree $T(\theta)$, let $R_l$ denotes its $l$-th terminal leaf. Since the space $\mathcal{X}$ is split into disjoint leaves by $T(\theta)$, we know for any $x \in \mathcal{X}$, there is exactly one leaf containing $x$. We let the index of the leaf be $l(x;\theta)$ and we say $x \in R_{l(x;\theta)}$.
Then for any single tree $T(\theta)$, the prediction on any data point $x \in \mathcal{X}$ is $ \sum_{i=1}^n w(X_i, x; \theta) Y_i$ where
Then a random forest containing $m$ trees formulates a prediction of $\mathbb{E}[Y|X=x]$ as \[ \sum_{i=1}^n w(X_i, x) Y_i \] where
From now on, we call the weight $w(X_i,x)$ in (ref) random forest weight. One can easily show that $\sum_{i=1}^n w(X_i,x) = 1$. The above representation of the random forest prediction of the mean can be equivalently obtained as a solution to the following least-squares optimization problem \[ \min_{\lambda \in \mathbb{R}} \sum_{i=1}^n w(X_i,x) (Y_i - \lambda)^2. \] Therefore, a least-squares regression adjustment, as the above, is equivalent to breiman2001random representation of random forests. However, when we move to estimation quantities that are not the mean, the latter representation is very powerful. Namely, a quantile random forest of meinshausen2006quantile can be seen as a quantile regression adjustment li2017forest, i.e., as a solution to the following optimization problem \[ \min_{\lambda \in \mathbb{R}} \sum_{i=1}^n w(X_i, x) \rho_{\tau}(Y_i - \lambda), \] where $\rho_\tau$ is the $\tau$-th quantile loss function, defined as $\rho_{\tau}(u) = u(\tau - \operatorname*{\mathbbm{1}}(u < 0))$. Local linear regression adjustment was also recently utilized in athey2016generalized to obtain a smoother and more poweful random forest algorithm.
Let us consider the case of no censored observations. Full data serve as a motivation for developing suitable estimating equations. Following the regression adjustment reasoning, for the case of fully observed data, we could estimate the $\tau$-th quantile of $T_i$ at $x$, denoted as $q_{\tau, x}$, as a solution to \[ \min_{q \in \mathbb{R}} \sum_{i=1}^n w(X_i, x) \rho_{\tau}(T_i - q). \] Equivalently, such estimate would solve the following estimating equations
where the second equality is true because $\sum_{i=1}^n w(X_i, x) = 1$.
For simplicity and better illustration of the idea, we first assume the latent variable $T_i$ has the same conditional probability in a neighborhood $R_x$ of $x$. Out of the $n$ data points, assume $\{X_1, \cdots, X_k\} \subset R_x$ and $w(X_i, x) = 1/k$ when $X_i \in R_x$ and $0$ otherwise. Now the estimating equation becomes
Now conditional on $\{x\} \cup \{X_i\}_{i=1}^k$,
which will be zero at $q^*$ where $\mathbb{P}(T > q^* | x) = 1-\tau$, that is, when $q^* = q_{\tau,x}$
Let's now consider the case of right-censored setting, where we further have the censoring variable $C_i$, which is independent of $T_i$ conditional on $X_i$, and we could only observe $Y_i = \min\{T_i, C_i\}$ and censoring indicator $\delta_i = \operatorname*{\mathbbm{1}} (T_i \le C_i)$. In order to estimate $q_{\tau,x}$, we cannot simply replace $T_i$ with $Y_i$ in (ref) as the $\tau$-th quantile of $T_i$ is no longer the $\tau$-th quantile of $Y_i$ because of the censoring. However, we can observe and utilize the following relationship
where $G(u|x)$ is the survival function of $C_i$ at $x$. That is to say, the $\tau$-th quantile of $T_i$ is actually the $1 - (1-\tau)G(q_{\tau,x}|x)$-th quantile of $Y_i$ at $x$. Now, we define a new estimating equation that resembles (ref) as follows
If we substitute $G(q|x)$ with $G(q_{\tau,x}|x)$, an intuitive explanation for (ref) is that as the $\tau$-th quantile of $T_i$ happens to be the $1 - (1-\tau)G(q_{\tau,x}|x)$-th quantile of $Y_i$ at $x$, instead of estimating the former which is not available because of the censoring, we turn to estimate the later one. Namely, the conditional expectation, $\mathbb{E}[S^o_k(q)]$, will still be zero at the same root $q^*$ for (ref).
The survival function $G(\cdot|x)$ can be estimated by a consistent estimate, for example the Kaplan-Meier estimator $\hat{G}(\cdot|x)$ using $\{Y_i\}_{i=1}^k$ and $\{\delta_i\}_{i=1}^k$, and we can then define
In the previous subsection, we made an assumption that $\mathbb{P}(T|X) = \mathbb{P}(T|x)$ for all $X \in R_x$, where $R_x$ is a neighborhood of $x$. But in reality, this assumption is not always true, and that is why $w(X_i, x)$ plays an important rule in our final estimator, as it “corrects" the empirical probability of each $T_i$ at $x$.
For example, say we have $n$ data points $\{(X_i, T_i)\}_{i=1}^n$ and have two cases: (1) at all $X_i$'s we have the same conditional probability of $T$, i.e. $\mathbb{P}(T|X_i) = \mathbb{P}(T|X_j)$ for all $i,j$; (2) $T$ has different conditional probabilities at different locations. In the setting (1), $X_i$'s become irrelevant and the point mass on each $T_i$ is $1/n$. We share the mass uniformly to the $n$ points $T_i$'s as they are equally important. When $n \to \infty$, it is known that for any $q$,
However, in the case (2), the convergence (ref) is no longer valid. We cannot simply put a mass $1/n$ on each $T_i$ because the probability of $T_i$ showing up at $X_i$ could be severely different than the probability it shows up at $x$. An extreme example is when $\mathbb{P}(T|x) = \text{Unif}(x-1,x+1)$. Then if $|X_i - x| > 1$, any $T_i$ showing up at $X_i$ should not even be counted when estimating $\mathbb{P}(T|x)$ because $\mathbb{P}(T_i|x) = 0$. In another word, we should give $T_i$ mass $0$ instead of $1/n$.
Therefore, a measure of “similarity" between points $X_i$ and $x$ needs to come into play, because we can no longer uniformly distribute the mass; observe that some $T_i$'s are more important than others for estimating $\mathbb{P}(T|x)$. For instance, if $X_i = x + 0.01$ and $X_j = x + 2$ in the previous example, then $T_i$ should be assigned much more weight than $T_j$.
Now let $w(X_i,x)$ denote the weight (mass) we assign to $T_i$ when we are estimating $\mathbb{P}(T|x)$. In the setting (1), we just have $w(X_i,x) = 1/n$ uniformly. But in the setting (2), we should have $w(X_i, x) > w(X_j, x)$ when $X_i$ is more similar to $x$ than $X_j$ in some sense. Therefore, the estimator for $\mathbb{P}(T \le q|x)$ is then
and it becomes clear that a proper weight $w(X_i,x)$ needs to satisfy:
One may naively think that any fixed Kernel weights, $K(X_i, x)$, could be a suitable choice. However, they would not be able to satisfy the second condition in (ref) for any distribution $\mathbb{P}(T | x)$. Fortunately, as shown in meinshausen2006quantile, the data-adaptive random forest weight $w(X_i, x)$ introduced in Section (ref) perfectly satisfy both conditions in (ref). And therefore going back to (ref), we now consider estimating equations,
when $n \to \infty$. Then following the same logic of how we get (ref), a heuristic extension of (ref) to censoring case will be
In the simplified example in Section (ref), we assume that $Y$ has the same conditional probability $\mathbb{P}(Y|X)$ in a neighborhood $R_x$ of $x$, and hence, we can estimate $G(q|x)$ by Kaplan-Meier estimator kaplan1958nonparametric (assuming no tied events)
where $k = |R_x|$. In the more complex case like in Section (ref), many consistent estimators for the conditional survival functions exist. For example, the nonparametric estimator proposed by beran1981nonparametric
is shown to be consistent beran1981nonparametric,dabrowska1987non,dabrowska1989uniform,gonzalez1994asymptotic,akritas1994nearest,li1995approach,van1996uniform. Here, $W_i(x,a_n)$ are the Nadaraya-Watson weights
$K(\cdot)$ is a known kernel and $\{a_n\}$ is a bandwidth sequence tending to zero as n tends to infinity. We can then simply use $\tilde{G}(q|x)$ as $\hat{G}(q|x)$ in (ref).
However, since we already have an adaptive version of kernel -- the random forest weights $w(X_i, x)$, we will propose the following two new estimators for $G(q|x)$.
\paragraph{Kaplan-Meier using nearest neighbors.} The first estimator is resembles that of (ref). We first find the $k$ nearest neighbors of $x$ according to the weights $w(X_i, x)$. Denoting these points as a set $N_x$, then we can simply use the Kaplan-Meier estimator on $N_x$
Here, the number of neighbors $k$ will be a tuning parameter.
\paragraph{Beran estimator with random forest weights.} In the second proposal, we will replace the Nadaraya-Watson weights in (ref) with random forest weights and get
One could observe that (ref) is a special case of (ref) when the weight $w(X_i, x) = 1/k$ for $X_i \in R_x$ and $0$ otherwise.
Finally, we summarize our main algorithm in Algorithm (ref). The details for choosing the candidate set $\mathcal{C}$ is in Section (ref). The choice to minimize the absolute value of $S_n(q;\tau)$ is arbitrary. The goal is to find the approximate root of $S_n(q;\tau) = 0$.
In this section, we will assume the random forest has terminal node size $m$, feature vector $X_i \in \mathbb{R}^{p}$, sample size is $n$, and $k$ nearest neighbors are chosen if using (ref).
The step (ref) in Algorithm (ref) involves of finding the $q^*$ in a candidate set $\mathcal{C}$ that sets the estimating equation $S_n(q; \tau)$ closest to zero. We simply evaluate the function $S_n(q;\tau)$ for all possible $q$ in $\mathcal{C}$ and find the minimum point. Note that for any fixed $\tau$, $S_n(q; \tau)$ is a step function in $q$ with jumps at $Y_i$'s because the discontinuities only happen at $Y_i$'s for $\hat{G}(q|x)$ (both (ref) and (ref)) and $\sum_{i=1}^n w(X_i, x) \operatorname*{\mathbbm{1}}(Y_i > q)$. Therefore, the candidate set $\mathcal{C} \subset \{Y_i\}_{i=1}^n$, and $|\mathcal{C}| = n$ in the worst case.
But in fact, for any fixed $x$, only $Y_i$'s with the corresponding feature vector $X_i \in R_x$ (ref) or with $w(X_i, x) > 0$ (ref) will be jump points, and hence, we can refine $\mathcal{C} = \{Y_i: X_i \in R_x\}$ for (ref) or $\mathcal{C} = \{Y_i: w(X_i, x) > 0\}$ for (ref). We then have the following theorem. The proof is given in the Appendix (ref).
In this section, we will show that for any fixed $\tau \in (0,1)$, $S_n(q; \tau)$ in (ref) will converge in probability to $(1-\tau)G(q|x) - \mathbb{P}(Y_i > q)$ uniformly for $q$.
We note that Condition (ref) is a very primitive condition on the distribution of the covariates. It is satisfied for example for Gaussian distribution and more broadly for most symmetric, continuous distributions with unbounded support. The case of bounded or discrete covariates is beyond the scope of the current work.
The two requirements of Condition (ref) are also required in meinshausen2006quantile (see Assumptions 2 and 3 therein). This condition states that the leaf node size of each tree should increase with the sample size $n$, but at a slower rate. Intuitively, first, the trees that we are using need to be shallow (i.e., with large leaves) in order to estimate a more complex model, reliably. Secondly, there can not be leaves with no samples, i.e., each leaf must be large enough to capture the local estimating equations more adequately. Our experiments also justify the necessity of Condition (ref), as the performance of our model, will deteriorate if we keep a small leaf node size but increase the sample size. We will talk about this in detail in Section (ref).
We note that Condition (ref) appears in all existing work related to quantile regression and inference thereafter.
Conditional independence of $T$ and $C$ is a very standard assumption and can be traced back to robins1992semiparametric among other works.
Condition (ref) is satisfied, for example, by the Kaplan-Meier estimator (ref) dabrowska1989uniform. Please take a look at Figure (ref) and Figure (ref) where we compare finite sample properties of the newly introduced estimators (ref) and (ref). We observe that the new distributional estimators are more adaptive and yet seemingly inherit consistency to that of the traditional KM estimator.
We proceed to showcase asymptotic properties of the proposed estimating equations. We begin by illustrating a concentration of measure phenomenon for the introduced score equations.
Next, we present our main result that illustrates an asymptotic consistency of the proposed conditional quantile estimator. The proof is given in Appendix (ref).
In this section, we will compare our model, censored forest regression (crf) with generalized random forest (grf) athey2016generalized, quantile random forest (qrf) meinshausen2006quantile and random survival forest (rsf) hothorn2005survival on simulated and real data sets.
On the simulated data sets, we will apply qrf and grf to the censored data directly, and get biased models which we denote by qrf and grf, respectively. We also apply qrf and grf to the data with uncensored responses, and call the resulted models \textit{qrf-oracle} and \textit{grf-oracle}.
Throughout this section, we fix the number of trees for each forest to be 1000. The only tuning parameter we have is the node size of each tree. All other parameters are kept as default.
In this section, we generate the true response, i.e., latent, variables $T_i \sim \textrm{Unif}(0,1)$, and consider the censoring variables $C_i \sim \mathcal{N}(0.8, 0.2^2)$. The censored responses, i.e. observed responses, are then taken as $Y_i = \min(T_i, C_i)$. We compare the estimating equation on the latent variables $T_i$
to the estimating equation of our proposed algorithm
where $\hat{G}(q)$ is the one-dimensional Kaplan-Meier estimator for the survival function of censoring variable $C$. The results are shown in Figure (ref). We consider $\tau=0.5$ but the results persist for many other choices of $\tau$.
In Figure (ref) we present the two estimating equations as functions of $q$ and illustrate that the solutions to $U_1(q)=0$ and $U_2(q)=0$ are closer and closer together when the sample size grows. The solution for $U_1(q)=0$ can be treated as an oracle solution where the oracle observes “uncensored" (true) response variable. Figure (ref) therefore indicated that the root of our method's estimating equation is very close to the oracle root and that we are therefore finding a good approximation to the unknown parameter of interest.
In this section, we will compare different forest algorithms on simulated data sets including accelerated failure time model (AFT) and many non-parametric censored regression models.
We simulate data from an one-dimensional AFT model
where $X \sim \textrm{Unif}(0,2)$ and $\epsilon \sim \mathcal{N}(0, 0.3^2)$. Then the censoring variable $C \sim \textrm{Exp}(\lambda = 0.08)$, and the observed response $Y = \min(T, C)$ and the censoring indicator $\delta = \operatorname*{\mathbbm{1}}(T \le C)$. The average censoring rate is about $20\%$. The number of training data, validation data and test data are all 300. All the forests consist of 1000 trees. The node size of each forest is determined by cross--validation. We plot out one set of training data and the corresponding quantile predictions for $\tau = 0.3, 0.5, 0.7$ on a set of test data in Figure (ref). We only show the results of crf, grf, and grf-oracle because in one dimension, qrf's performance is visually indistinguishable from grf. There we observe a consistency of our method as well as superior behavior to the competing methods. Namely, the generalized random forest that ignores the censoring component of the data, incurs large bias; due to the right censoring, bias is larger for lower values of the quantiles. We observe that the proposed {\it crf} follows closely the oracle estimator and is extremely close to the true quantile regardless of the $\tau$ in the study.
Moreover, we proceed further and for a set of values $\tau \in \{0.1,0.3,0.5,0.7,0.9\}$, we repeat the process 40 times, and for each time, we calculate the MSE and MAD between the estimated quantiles and the true quantiles, and the $\tau$-th quantile loss. To be more specific, let $T_i$ be the response in test set (all uncensored), $Q^{\tau}_i$ be the true $\tau$-th quantile, and $\hat{Q}^{\tau}_i$ be the estimated quantile, then
The reason we use $L_{quantile}$ to measure the quality of quantile predictions is that, by meinshausen2006quantile, the $\tau$-th quantile of $T$ at $x$ equals to $\operatorname*{arg\,min}_{q \in \mathbb{R}} \mathbb{E}[\rho_{\tau}(T - q) | X=x]$. The results are illustrated in Figure (ref) where besides the abose three measures we compare the concordance index (C-index) harrell1982evaluating, which is related to the area under the ROC curve heagerty2005survival. It estimates the probability that, in a randomly selected pair of cases, the case that fails first had a worse predicted outcome. In ishwaran2008random, they use the ensemble mortality as the predictive outcome for their random survival forest, and the predicted survival time for random forest regression. For our method crf and the other two methods, qrf and grf, we will use the $\tau$-th conditional quantile as the predicted outcome. Since the outcomes will be different for different $\tau$, we report the results for all $\tau \in \{0.1,0.3,0.5,0.7,0.9\}$.
In Figure (ref) we observe an oracle like behavior of the proposed {\it crf} method in terms all four measures of the quality of estimation and/or prediction. Namely, we observe that MAD, MSE and quantile losses are extremely small whereas C-index is high and all are close to the corresponding oracle estimators (colored purple and blue). Moreover, we observe that the proposed {\it crf} method, although not primarily build for the hazard rates, is even better than survival random forest: see for example discrepancies between red and brown boxplots in the last row of Figure (ref) where the larger the C-index is the better the method is.
In this section, we will compare the two different conditional survival function estimators (ref) and (ref). We generate training data and test data from the one-dimensional AFT model defined in the previous section, but with two different censoring rate:
We then choose four test points $\{x_1=0.4, x_2=0.8, x_3=1.2, x_4=1.6\}$, and then plot out the conditional survival function estimators $\hat{G}(q|x_i)$ by the two different methods (ref) and (ref) on these four points. The results are shown in Figure (ref) and (ref) for three different training sample sizes $n \in \{300, 2000, 5000\}$. For the nearest neighbor estimator (ref), we set the number of neighbors to be $n/10$, which is also the node size we choose.
We can observe that when $n$ increases, two curves become closer and are both good approximations of the true survival curve. But the first method (ref) does have an extra tuning parameter $k$ -- the number of nearest neighbors, so in the experiments, we always choose to use the second estimator (ref), which is more adaptive and parameter free.
Note that the estimated survival function will degenerate at the tail of the distribution when the test point $x$ is small; see the first two columns in Figure (ref) and (ref). This is a common phenomenon even for the regular KM estimator because there is no censored observations beyond some time point. In the AFT model, the conditional distribution of the latent variable depends on the location $x$. When $x$ is small, the conditional mean of $T$ is also small, and we could not observe most of the censoring values where $C_i > T_i$, leading to degenerated survival curves. However, if we continue increasing the sample size $n$, we should be able to recover the entire curve even for smaller $x$. In fact, when we increase the censoring level from $20\%$ (Figure (ref)) to $50\%$ (Figure (ref)), we find that both estimators give better performance because we can observe more censored values.
Since our forest regression method crf is nonparametric and does not rely on any parametric assumption between response and explanatory variables, it can be used to estimate quantiles for any general model $T = f(X) + \epsilon$. In this section, we let $f(x) = \sin(x)$ and have the model
where $X \sim \textrm{Unif}(0,2 \pi)$ and $\epsilon \sim \mathcal{N}(0, 0.3^2)$. Then the censoring variable $C \sim 1 + \sin(X) + \textrm{Exp}(\lambda = 0.2)$, and the responses $Y = \min(T, C)$. All the other settings are the same as in Section (ref). We plot out the training data and the quantile predictions for $\tau = 0.3, 0.5, 0.7$ in Figure (ref). The censoring level is about $25\%$. We observe that for all $\tau \in \{0.3,0.5,0.7\}$, crf can produce comparable quantile predictions to grf-oracle. Especially when $\tau=0.3$, the quantile prediction by $grf$ (blue dotted curve) severely deviates from the true quantile, while our method crf can still predict the correct quantile and performs as good as grf-oracle. We want to emphasize that grf-oracle uses the latent responses $T_i$ while our method only uses the observed responses $Y_i$ and censoring indicators $\delta_i$. We then repeat the experiments for 40 times and report the box plots in Figure (ref). Again we can see that for all quantiles, our method crf behaves almost as good as \textit{qrf-oracle} and \textit{grf-oracle}, and consistently better than \textit{qrf} and \textit{grf}. For example the order of magnitude of our error is twice and sometimes more than two times smaller than that of quantile or generalized random forest.
In this section, we test our algorithm on a multi-dimensional AFT model
where $\bm{\beta} = (0.1,0.2,0.3,0.4,0.5)$, $X_{\cdot,j} \sim \textrm{Unif}(0, 2)$, and $\epsilon \sim \mathcal{N}(0, 0.3^2)$. The censoring variable $C \sim \textrm{Exp}(\lambda = 0.05)$, and $Y = \min(T, C)$. The censoring level is about $22\%$. The number of training data is 500 and the number of test points is 300. All the forests consist of 1000 trees. The result is in Figure (ref). Our model crf still outperforms qrf and grf significantly, and is comparable to qrf-oracle and grf-oracle.
In this section, we construct a complex model
where $X_{\cdot,j}$ stands for $j$-th dimension of $\bm{X} \in \mathbb{R}^5$, and $\epsilon \sim \mathcal{N}(0,0.3^2)$. Then we consider a censoring variable independent of $\bm{X}$ and $T$: $C \sim \textrm{Exp}(\lambda = 0.015)$. The result is summarized in Figure (ref). Although the model is highly non-linear our method is able to capture that and estimate the quantile of interest consistently. However, other methods are unable to be consistent. Behavior matches that of linear models closely. In this case we notice that the oracular generalized random forests have much better performance; note that since we are estimating equations, our algorithm can be easily tweaked to match {\it grf-oracle} -- the only change would be to design the splits in the initial tree construction to match our estimating equations, use subsampling and save a separate sample for minimizing the equations themselves.
In this section, we investigate the impact of the choice of the node size on different methods. The data we use will be generated from the one-dimensional and multi-dimensional AFT and Sine models as defined in the previous sections. We increase the node size from 5 to 60 with step size of 5.
\paragraph{One-dimensional AFT and Sine models} The result of sine model is summarized in Figure (ref). One can see that for both qrf and our model, crf, the quantile loss will first decrease when node size increases. It attains minimum around node size of 30. However, for grf, its quantile loss is almost monotonically increasing, and attains minimum at node size of 5. But both qrf-oracle and grf-oracle can attain the best quantile loss of about 0.125. And one impressive observation is that our model, crf, almost performs the same as \textit{qrf-oracle} for all node sizes. Similar conclusion can be made from the AFT result which is in Figure (ref).
\paragraph{Multi-dimensional AFT model} From Figure (ref), we observe that qrf, qrf-oracle and grf-oracle all give similar results. The performance of our model crf is only slightly worse than qrf-oracle, but is even better than grf-oracle.
\paragraph{Multi-dimensional complex model} The result is summarized in Figure (ref). The censoring level is about $25\%$. From the figure, we observe that the behavior of crf is still only slightly worse than qrf-oracle. In this experiment, grf-oracle behaves the best. All of crf, qrf-oracle and grf-oracle are significantly better than the biased models \textit{qrf} and \textit{grf}. When $\tau=0.7$, \textit{grf} behaves slightly better than \textit{qrf-oracle} when node size is small. The reason is that the conditional quantiles of $Y$ and $T$ are closer when $\tau$ is larger, and \textit{grf} is more stable and smooth on the data in this experiment. But we still observe that the performance of \textit{crf} and \textit{qrf-oracle} are very close.
In this section, we compare our censored forest (crf) with quantile random forest (qrf) meinshausen2006quantile and generalized forest (grf) athey2016generalized on real datasets. In order to evaluate the performances unbiasedly, we manually add censoring to the data. In addition, we apply qrf and grf to the data without censoring and we call the resulted models qrf-oracle and grf-oracle, respectively.
For all these methods, bagged versions of the training data are used for each of the $1000$ trees. We use 5-fold cross validation to select the best node sizes for different methods. For all the other parameters, we keep the default settings.
\paragraph{Datasets} We use datasets BostonHousing, Ozone from the R packages mlbench and alr3. For all the datasets, we sample censoring variables from Exponential distributions with $\lambda$ set so that the censoring level is roughly $20\%$. For BostonHousing dataset, we set $\lambda = 0.01$. For Ozone, $\lambda = 0.025$. For \textit{Abalone} dataset, we random sample $1000$ observations and take the log-transformation of the response variable \textit{rings}. We then set $\lambda = 0.10$.
\paragraph{Evaluation} For each dataset, we train our model on bootstrapped version of the data, and test the performance on out-of-bag observations. This process is repeated for 40 times, and we calculate the mean and standard deviation of the prediction errors. In our context, the prediction error is measured by the $\tau$-th quantile loss for $\tau$-th quantile estimation. The results are illustrated in Figure (ref).
On all data sets, our proposed method behaves better than quantile forest and generalized forest in terms of quantile losses. Especially when $\tau = 0.1$, $0.3$ or $0.5$, the performance of our method is significantly better than qrf and grf, and is even comparable to that of oracle qrf and grf. It agrees with our observation in the one-dimensional example (Figure (ref) and (ref)). While estimating larger quantiles, the true $\tau$-th quantile of $T_i$ and $Y_i$ are close, and hence the performance of all five models are similar. But when $\tau$ is small, the $\tau$-th quantile of $T_i$ and $Y_i$ are different because of the censoring, and in this case, our model has superior advantage and find the true quantiles of $T_i$ almost as good as the oracle methods.
\paragraph{Nodesize} For each dataset, we train different models using different nodesizes and compare the performance. For each node size, we bootstrap the data and repeat the experiments for 20 times, and we calculate the mean and standard deviation of the quantile predictions for quantiles $\tau = 0.3$, $0.5$, and $0.7$. The result is in Figure (ref). We observe that our method, crf, is uniformly better than qrf and grf, proving that crf is able to correct the bias introduced by censoring. Moreover, the quantile loss of crf is always competitive to that of qrf-oracle and \textit{grf-oracle}, and is remarkably always better than \textit{grf-oracle}, only slightly worse than \textit{qrf-oracle}.
All the forest methods can be used to get $95\%$ prediction intervals by predicting the $0.025$ and $0.975$ quantiles of the true response variable. Then for any location $x \in \mathcal{X}$, a straightforward confidence interval will be $[Q(x;0.025), Q(x;0.975)]$. The result is illustrated in Figure (ref) for the case of univariate censored sine model. For each data set, we bootstrap the data and calculate the $0.025$ and $0.975$ quantile for the out of bag points. Then for each node size, we repeat this process for 20 times and calculate the average coverage rate of the confidence intervals.
We observe that in all of the cases, our method crf and qrf-oracle give the coverage closest to $95\%$. As can be seen from Figure (ref), both qrf and grf perform much worse on predicting lower quantiles. They tend to under-estimate the lower quantiles and hence make the confidence intervals much wider than the true ones. Nevertheless, as seen in Figure (ref), out method has great coverage that is never below $95\%$, indicating certain efficiency of the proposed method.
In this article, we introduced a censored quantile random forest, a novel non-parametric method for quantile regression problems that is integrated with the censored nature of the observations. While preserving information carried by the censored observations, the novel estimating equations maintain the flexibility of general random forest approaches. The estimating equations are equipped with both censoring information as well as random forest information simultaneously.
The statistical mechanism, particularly asymptotic normality, of ensemble tree-based methods is still not fully understood. Some insightful discussions and first steps towards this goal can be found in zhu2012recursively. As noted in athey2016generalized even consistency properties of censored based methods present significant theoretical challenges. Generalization of the results of athey2016generalized to the case of censored observations involves non-trivial extensions of the theoretical advancements introduced therein. A generalization would require not only adaptation to censoring but as well as to an infinite-dimensional nuisance parameter; current results focus on finite dimensional nuisance parameters. Once the latter generalization is achieved then, utilizing our estimating equations, will allow for the generalization of the former.
One of the promising applications of the introduced method is in the estimation of heterogeneous treatment effects when the response variable is censored; treatment discovery with right-censored observations is important and yet poorly understood research area. Equipping this literature with a fully non-parametric approach would lead to a significant broadening of the now more known parametric approaches. Our method appears to be a nice fit for this setting; observe that our estimating equations can be easily replaced with another kind that targets treatment effects directly.