Survival ensembles

Survival ensembles TORSTEN HOTHORN∗ ... Seminar fur Statistik, ... ally efficient form of regularization applicable to a wide class of linear models in...

0 downloads 125 Views 673KB Size
Biostatistics (2006), 7, 3, pp. 355–373 doi:10.1093/biostatistics/kxj011 Advance Access publication on December 12, 2005

Survival ensembles TORSTEN HOTHORN∗

¨ PETER BUHLMANN Seminar f¨ur Statistik, ETH Z¨urich, CH-8032 Z¨urich, Switzerland SANDRINE DUDOIT Division of Biostatistics, University of California, Berkeley, 140 Earl Warren Hall, #7360, Berkeley, CA 94720-7360, USA ANNETTE MOLINARO Division of Biostatistics, Epidemiology and Public Health, Yale University School of Medicine, 206 LEPH, 60 College Street, PO Box 208034, New Haven CT 06520-8034, USA MARK J. VAN DER LAAN Division of Biostatistics, University of California, Berkeley, 140 Earl Warren Hall, #7360, Berkeley, CA 94720-7360, USA

S UMMARY We propose a unified and flexible framework for ensemble learning in the presence of censoring. For right-censored data, we introduce a random forest algorithm and a generic gradient boosting algorithm for the construction of prognostic and diagnostic models. The methodology is utilized for predicting the survival time of patients suffering from acute myeloid leukemia based on clinical and genetic covariates. Furthermore, we compare the diagnostic capabilities of the proposed censored data random forest and boosting methods, applied to the recurrence-free survival time of node-positive breast cancer patients, with previously published findings. Keywords: Censoring; Cross-validation; Ensemble methods; IPC weights; Loss function; Prediction; Prognostic factors; Survival analysis.

∗ To whom correspondence should be addressed. c The Author 2005. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: [email protected] 

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

Institut f¨ur Medizininformatik, Biometrie und Epidemiologie, Friedrich-Alexander-Universit¨at Erlangen-N¨urnberg, Waldstraße 6, D-91054 Erlangen, Germany [email protected]

356

T. H OTHORN ET AL . 1. I NTRODUCTION

In survival time studies, models regressing the time to event on a set of covariates, i.e. variables expected to be associated with the disease the patient suffers from, are the basis of prognostic and diagnostic modeling. The specification and estimation of such models are complicated by the fact that often only incomplete information about the response variable is available due to censoring. The most widely used representative of regression methods for censored data is the Cox proportional hazards model (Cox, 1972), which addresses the censoring problem by maximizing the partial likelihood while leaving the baseline hazard unspecified under the proportional hazards assumption. In order to motivate the methodology proposed in this paper, it is helpful to classify existing approaches as addressing one of the following four problems.

Model assumptions. Many authors proposed flexible alternatives to the Cox model without assuming proportional hazards, such as (partially) nonlinear accelerated failure time models (Stute, 1999; Orbe et al., 2003), spline-based extensions (Gray, 1992; Kooperberg et al., 1996; LeBlanc and Crowley, 1999), fractional polynomials (Sauerbrei and Royston, 1999), and neural networks (Ripley et al., 2004). Dimensionality. Current research efforts have focused on data analysis problems with high-dimensional covariate spaces, mainly driven by the requirements of biological applications such as microarray gene expression profiling. In high-dimensional situations, Hastie and Tibshirani (2004) suggest a computationally efficient form of regularization applicable to a wide class of linear models including the Cox model, while Huang and Harrington (2005) investigate iterative partial least-squares fitting in accelerated failure time models. In contrast, dimension reduction techniques are studied by Li and Li (2004) and Bair and Tibshirani (2004), who advocate the application of low-dimensional compound covariates obtained from an unsupervised clustering of the covariates. Model selection and evaluation. Another important research problem concerns model selection and evaluation. While classical techniques like residual analysis (e.g. Therneau and Grambsch, 2000) and the detection of influential observations (Bedrick et al., 2002) have been translated into the context of survival analysis, specialized goodness of prediction measures, such as the Brier score for censored data (Graf et al., 1999), are a matter of debate (Henderson, 1995; Altman and Royston, 2000; Schemper, 2003). Although censoring induces nontrivial problems for the comparison of observed and predicted responses, such measures are important for cross-validation and other resampling-based model selection and evaluation techniques (Sauerbrei, 1999; van der Laan and Dudoit, 2003; Dudoit and van der Laan, 2005; Hothorn et al., 2005b). In this paper, we address the four aforementioned problems simultaneously, by applying the general estimation framework described in van der Laan and Robins (2003) in order to generalize ensemble learning techniques to censored data problems. The framework allows for the specification of regression models under complete information (full data world) for arbitrary loss functions. For estimation of the models under incomplete information (observed data world), a special weighting scheme ensures that observations likely to be censored are up-weighted compared to observations on patients likely to experience an event. As a consequence, in the absence of censoring, the models reduce to their counterparts known from the uncensored situation. Most importantly, the goodness of prediction of such models is easily evaluated

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

Connection between censored and uncensored methods. The establishment of a close connection between regression models for uncensored continuous response variables and models designed for censored data was motivated by the problem that the Cox model does not reduce to an ordinary linear regression model in the absence of censoring. Accelerated failure time models (e.g. James, 1998), such as the Buckley–James model (Buckley and James, 1979), do have this desirable property.

Survival ensembles

357

2. M ODEL The estimation problems to be solved are first defined in the full data world and are then mapped into the observed data world, i.e. in the presence of censoring, following van der Laan and Robins (2003). 2.1

Full data world

In an ideal world, we are able to observe random variables Z = (Y, X) with distribution function FY,X , where Y = log(T ) and T ∈ R+ denotes a survival time. The p-dimensional covariate vector X = (X 1 , . . . , X p ) is taken from a sample space X = X1 ×· · ·×X p . We assume that the conditional distribution FY |X = FY | f (X) of the response Y given the covariates X depends on X through a real-valued function f : X → R. The regression function f , our parameter of interest, is an element of some parameter space  and has minimal risk   EY,X L(Y, f (X)) = L(Y, f (X)) dFY,X = min L(Y, ψ(X)) dFY,X ψ∈

for a suitable full data loss function L: R × R → R+ . Our principle aim is to estimate the regression function f . Usually, an estimate fˆ of f is computed via constrained minimization of the empirical risk defined by the full data loss function L. However, this minimization problem can only be solved when all quantities are observed. Naturally, this is not the case in the presence of censoring.

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

using cross-validation techniques based on well-known loss functions such as the quadratic or absolute loss (Keles¸ et al., 2004; van der Laan and Dudoit, 2003). The general estimation framework has recently been applied to longitudinal marginal structural models (Bryan et al., 2004), survival trees (Molinaro et al., 2004), and other estimation problems (Sinisi and van der Laan, 2004; van der Laan et al., 2004). Ensemble methods like bagging, random forest, and boosting yield flexible predictors for nominal and continuous responses and are known to remain stable in high-dimensional settings. Briefly, bagging and random forest average over predictions of simple models, so called base learners, which have been fitted to slightly perturbed instances of the original observations, whereas boosting is a functional gradient descent algorithm with iterative fitting of appropriately defined residuals. For a general overview, we refer to B¨uhlmann (2004) and references therein. Here, we extend the area of application of ensemble methods to survival analysis. We incorporate weights into random forest-like algorithms and extend gradient boosting in order to minimize a weighted form of the empirical risk. The published attempts to use ensemble techniques for modeling censored data are rather limited due to the difficulties induced by censoring. Ridgeway (1999) proposed a boosting algorithm minimizing the partial likelihood and Benner (2002) derived a boosting algorithm from the Brier score for censored data. A special aggregation scheme for bagging survival trees was studied by Hothorn et al. (2004). Ishwaran et al. (2004) constructed an random forest predictor of mortality from heart disease by averaging relative risk trees. Breiman (2002) introduced a software implementation of an random forest variant for censored data, but without a formal description of the methodology being available. Following the road map of van der Laan and Robins (2003) and van der Laan and Dudoit (2003), Section 2 defines the regression models and the corresponding risk optimization problems in the ‘full data world’ and sketches the general estimating framework in the ‘observed data world’. In Section 3, we propose both an random forest and a boosting algorithm for censored data. The advantages of our approaches are studied in Section 4 with respect to the stability and flexibility of predictions for patients suffering from acute myeloid leukemia (AML), based on high-dimensional covariates from gene expression profiling experiments and clinical data. Moreover, we focus on the diagnostic capabilities of flexible ensemble methods for data from the node-positive breast cancer patients.

358

T. H OTHORN ET AL . 2.2 Observed data world

A particular example for the nuisance parameter η will be given in Section 2.3. The basic idea is to minimize the empirical counterpart of EY˜ ,,X L(Y˜ , ψ(X)|η) with respect to the candidate estimators ψ ∈ , which is possible even in the imperfect observed data world. 2.3

Inverse probability of censoring weights

One approach for defining the observed data loss function L(Y˜ , ψ(X)|η) is the application of inverse probability of censoring (IPC) weights (van der Laan and Robins, 2003), for which the nuisance parameter η is given by the conditional censoring survivor function G: L(Y˜ , ψ(X)|G) = L(Y˜ , ψ(X))

 . G(T˜ |X)

Basically, the full data loss function is weighted by the inverse probability of being censored after time T˜ given the covariates X. The inverse probability G(T˜ |X)−1 exists because G(T˜ |X)  G(T |X) > 0 by assumption. The corresponding empirical risk is the weighted average Eˆ Y˜ ,,X L(Y˜ , ψ(X)|G) = n −1

n  i=1

ˆ = n −1 L(Y˜i , ψ(Xi )|G)

n  i=1

L(Y˜i , ψ(Xi ))

i ˆ T˜i |Xi ) G(

(2.1)

and the regression function estimator fˆ is derived by (constrained) minimization of (2.1) with respect to the candidate estimators ψ ∈ . Note that the conditional censoring survivor function G is typically ˆ A Kaplan–Meier estimate Gˆ is the simplest choice unknown and needs to be replaced by an estimate G. but other procedures, based for example on a Cox model, are appropriate. For convenience, let w = ˆ T˜i |Xi )−1 , denote the IPC weights. Other choices of the observed data (w1 , . . . , wn ), where wi = i G( loss function are possible, such as that based on doubly robust inverse probability of censoring (DR-IPC) weights (van der Laan and Robins, 2003; Molinaro et al., 2004).

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

In realistic set-ups we only observe random variables O = (Y˜ , , X), where Y˜ = log(T˜ ) for time to event T˜ = min(T, C) and censoring indicator  = I (T  C), from some distribution FY˜ ,,X . We assume that the conditional censoring distribution P(C  c|Z) only depends on the covariates, that is P(C  c|Z) = P(C  c|X), or, equivalently, that survival time T and censoring time C are conditionally independent given the covariates X. This assumption implies the ‘coarsening at random’ (CAR) assumption on the censoring mechanism (for details, we refer to van der Laan and Robins, 2003, Section 1.2.3). Furthermore, for the corresponding conditional censoring survivor function G(c|X) = P(C > c|X), we assume that G(T |X) is strictly greater than zero almost everywhere with respect to the full data distribution FY,X . The implications of this assumption are discussed in Section 5. The parameter space  is the function space of all candidate estimators ψ: X → R for the regression function f . For an observed learning sample of n independent and identically distributed observations L = {Oi = (Y˜i , i , Xi ); i = 1, . . . , n}, we cannot evaluate the full data loss function L(Y, ψ(X)) for the censored patients. Consequently, we cannot minimize the corresponding empirical risk defined in terms of the full data loss function L(Y, ψ(X)) directly. The methodology presented in van der Laan and Robins (2003) and van der Laan and Dudoit (2003) solves this problem by replacing the full data loss function L(Y, ψ(X)) by an observed data loss function L(Y˜ , ψ(X)|η) with nuisance parameter η, where the risks of both the loss functions coincide for all candidate estimators ψ ∈ :   EY,X L(Y, ψ(X)) = L(Y, ψ(X)) dFY,X = L(Y˜ , ψ(X)|η) dFY˜ ,,X = EY˜ ,,X L(Y˜ , ψ(X)|η).

Survival ensembles

359

3. E NSEMBLE LEARNING We present two algorithms pursuing some regularized minimization of (2.1): random forest and gradient boosting for censored data. The random forest approach seeks to minimize the empirical risk indirectly via a stabilization of randomized base learners fitted on perturbed instances of the learning sample L. In contrast, gradient boosting employs a functional gradient descent algorithm for minimizing the empirical risk (2.1). 3.1

Random forest

Algorithm: Random Forest for Censored Data Step 1 (Initialization). Set m = 1 and fix M > 1. Step 2 (Bootstrap). Draw a random vector of case counts vm = (v m1 , . . . , v mn ) from the multinomial  n −1 distribution with parameters n and w. i=1 wi Step 3 (Base Learner). Construct a partition πm = (Rm1 , . . . , Rm K (m) ) of the sample space X into K (m) cells via a regression tree. The tree is built using the learning sample L with case counts vm , i.e. is based on a perturbation of the learning sample L with observation i occurring v mi times. Computational details are given below. Step 4 (Iteration). Increase m by one and repeat Steps 2 and 3 until m = M. Prognostic modeling is our main concern, i.e. we are interested in estimating the (log)-survival time fˆ(x) for a patient with covariate status x. The predicted status of the response variable is computed based on prediction weights ai (x) =

M  m=1

v mi

K (m)

I (Xi ∈ Rmk and x ∈ Rmk );

i = 1, . . . , n.

k=1

The prediction weight ai (x) measures the ‘similarity’ of x to Xi (i = 1, . . . , n) by counting how many times the value x falls into the same cell as the ith observation in the learning sample. The prediction fˆ(x) can now be computed as the solution of Yˆ = fˆ(x) = argmin y∈R

n 

L(Y˜i , y)ai (x).

i=1

For quadratic loss L(Y, ψ(X)) = (Y − ψ(X))2 , the prediction is simply the weighted average of the observed (log)-survival times  n −1 n   ˆ |X = x) = fˆ(x) = Yˆ = E(Y ai (x) ai (x)Y˜i . i=1

i=1

The full data loss function can be evaluated here because, by definition, the weights wi , and thus the case counts v mi as well as the prediction weights ai (x), are zero for censored observations. The prediction weights approach is essentially an extension of the classical (unweighted) averaging of predictions extracted from each single partition (cf. Breiman, 1996) as used also in Hothorn et al. (2004).

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

From the observed learning sample L = {(Y˜i , i , Xi ); i = 1, . . . , n}, compute the weight vector w. Note that the learning sample can be thought to include the censored observations as well, but with wi = 0 iff i = 0. The random forest algorithm with weights w basically works by defining the resampling probability of observation i in terms of the corresponding weight wi .

360

T. H OTHORN ET AL .

In Step 3 of the algorithm, the partitions are usually induced by some form of recursive partitioning with additional randomization. This can be implemented by using only a small number of randomly selected covariates for further splitting of every node of the tree. Note that the proposed random forest for censored data reduces to the original random forest procedure (Amit and Geman, 1997; Breiman, 2001a) when all events have been observed. Conceptually, the algorithm is not restricted to (randomized) trees as base learners. Any other regression model can be applied as well. However, the prediction weights approach is only applicable when some form of recursive partitioning is used. For all other base learners, survival times need to be estimated via unweighted averages of the predictions similar to the original bagging approach. A practical drawback of the random forest algorithm for censored data is that out-ofbag predictions, and thus out-of-bag error rate estimates, cannot be computed when some observations are given a very large weight and are thus appearing in nearly every bootstrap sample.

In the full ndata world, the generic boosting algorithm sketched here can be applied to pursue minimization of i=1 L(Yi , ψ(Xi )) via functional gradient descent (for details, we refer to Friedman, 2001, and B¨uhlmann and Yu, 2003). Let U denote a pseudo-response variable. A base learner regressing the pseudoresponse U on the covariates X is denoted by h(·|ϑU,X ), where ϑU,X is a vector of parameters. Fitting the base learner can be performed by minimizing any loss function, for example solving the least-squares problem n  (Ui − h(Xi |ϑ))2 . (3.1) ϑˆ U,X = argmin ϑ

i=1

Algorithm: Generic Gradient Boosting for Uncensored Data Step 1 (Initialization). Define Ui = Yi (i = 1, . . . , n), set m = 0, and fˆ0 (·) = h(·|ϑˆ U,X ). Fix M > 1. Step 2 (Gradient). Compute the residuals

 ∂ L(Yi , ψ)  Ui = −  ˆ ∂ψ ψ= f m (Xi )

and fit the base learner h(·|ϑˆ U,X ) to the new ‘responses’ Ui as in (3.1). Step 3 (Update). Update fˆm+1 (·) = fˆm (·) + νh(·|ϑˆ U,X ) with step size 0 < ν  1, for example ν = 0.1. Step 4 (Iteration). Increase m by one and repeat Steps 2 and 3 until m = M. Note that, unlike the random forest algorithm, the number of iterations M is a tuning parameter which needs to be determined via cross-validation. Internal stopping criteria are available for special cases, which we will discuss in Section 3.4. 3.3

Gradient boosting—observed data world

In the observed data world, we cannot solve the least-squares problem (3.1) for fitting the base learner since we do not have access to Ui which is a function of the censored Yi . However, the right-hand side of (3.1) can be replaced by an empirical risk as in (2.1) and we then get the weighted least-squares problem  n  ∂ L(Y˜i , ψ)  2 ˜ ˜ ˆ wi (Ui − h(Xi |ϑ)) with pseudo-responses Ui = − . ϑU˜ ,X = argmin   ˆ ∂ψ ϑ i=1

Thus, the following algorithm can be applied to minimize the empirical risk (2.1).

ψ= f m (Xi )

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

3.2 Gradient boosting—full data world

Survival ensembles

361

Algorithm: Generic Gradient Boosting for Censored Data Step 1 (Initialization). Define U˜ i = Y˜i (i = 1, . . . , n), set m = 0, and fˆ0 (·) = h(·|ϑˆ U˜ ,X ). Fix M > 1. Step 2 (Gradient). Compute the residuals  ˜i , ψ)  ∂ L( Y U˜ i = −   ∂ψ

ψ= fˆm (Xi )

and fit the base learner h(·|ϑˆ U˜ ,X ) to the new ‘responses’ U˜ i by weighted least squares. Step 3 (Update). Update fˆm+1 (·) = fˆm (·) + νh(·|ϑˆ U˜ ,X ) with step size 0 < ν  1. The boosting estimator is fˆM (x) and the predicted (log)-survival time for an observation with covariate status x is Yˆ = fˆM (x). The algorithm proposed here reduces to the original form of gradient boosting in the absence of censoring. For quadratic loss L(Y, ψ(X)) = (Y − ψ(X))2 /2, the algorithm is obtained by residuals U˜ i = Y˜i − fˆm (Xi ) in the mth boosting iteration and we call this method L 2 -boosting for censored data. 3.4 Gradient boosting—choice of base learners and stop criterion The base learner h needs to be able to take weights w into account. Recursive partitioning procedures are popular choices of such base learners and the methodology of Molinaro et al. (2004) can be applied directly. B¨uhlmann and Yu (2003) suggested univariate smoothing splines: in each boosting iteration, one of the p covariates is selected, and the relationship between the residuals U and the selected covariate is modeled by a smoothing spline with low degrees of freedom. Another possibility, which is studied here, is the application of componentwise least squares (B¨uhlmann, 2006). This choice is computationally attractive and allows for the definition of an AIC-based internal stopping criterion. Let X( j) denote the design matrix associated with the jth covariate. When the jth covariate is a factor, the matrix X( j) is a dummy matrix. A column for the intercept term could be √ included. Let W denote the n × n diagonal matrix with diagonal elements Wii = wi , i = 1, . . . , n. Then, H( j) = X( j) ((WX( j) ) (WX( j) ))−1 (WX( j) ) W is the usual hat matrix of a simple linear model with covariate j alone. In the mth boosting iteration, we select the covariate with minimum empirical risk, i.e. km = argmin

n 

˜ i )2 , wi (U˜ i − (H( j) U)

j=1,..., p i=1

˜ = (U˜ 1 , . . . , U˜ n ) is the vector of pseudo-responses in the mth step. The fit in the mth step where U ˜ as introduced by can be written in terms of the boosting hat operator ( fˆm (X1 ), . . . , fˆm (Xn )) = Bm Y ˜ = (Y˜1 , . . . , Y˜n ) denotes the n-vector of responses extracted from B¨uhlmann and Yu (2003), where Y L. In the first boosting iteration, the boosting hat operator is B0 = νH(k0 ) and the update Step 3 can be written as Bm+1 = Bm + νH(km ) (In − Bm ), where the n × n matrix In denotes the identity matrix. This formulation of boosting in terms of a boosting operator opens up the way to an AIC-based internal stopping criterion (B¨uhlmann, 2006). The trace of the boosting operator Bm is interpreted as degrees of

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

Step 4 (Iteration). Increase m by one and repeat Steps 2 and 3 until m = M.

362

T. H OTHORN ET AL .

freedom. A corrected version of AIC can be computed as AIC(m) = log(σˆ 2 ) +

1 + trace(Bm )/n , 1 − (trace(Bm ) + 2)/n

where σˆ 2 = n −1

˜ i )2 , wi (Y˜i − (Bm Y)

i=1



where the weights have been rescaled to wi = wi ( boosting iterations is Mˆ = argminm=1,...,M AIC(m).

n 

i

wi )−1 n. An estimate of the optimal number of

4. I LLUSTRATIONS AND APPLICATIONS

4.1

Acute myeloid leukemia

The treatment of patients suffering from AML is determined by a tumor classification scheme taking the status of various cytogenetic aberrations into account. Bullinger et al. (2004) investigated an extended tumor classification scheme incorporating molecular subgroups of the disease obtained by gene expression profiling. A combination of unsupervised and supervised techniques is applied to define a binary outcome predictor (good versus poor prognosis) taking into account the expression measures of 133 selected genes which are represented by 149 complementary DNAs (cDNAs) . This binary surrogate variable is shown to discriminate between patients with short and longer survival in an independent sample of patients. Instead of using a binary variable summarizing expression levels of 149 cDNAs, random forest and L 2 -boosting are applied to construct predictors based on both the clinical data and the expression levels of the genes selected by Bullinger et al. (2004). The results reported here are based on clinical and gene expression data published online at http://www.ncbi.nlm.nih.gov/geo, accession number GSE425. The overall survival time and censoring indicator as well as the clinical variables age, sex, lactate dehydrogenase level, white blood cell count, and treatment group are taken from Supplementary Table 1 in Bullinger et al. (2004). In addition, this table provides two molecular markers, the fms-like tyrosine kinase 3 and the mixed-lineage leukemia gene, as well as cytogenetic information helpful to define a risk score (‘low’: karyotype t(8;21), t(15;17) and inv(16); ‘intermediate’: normal karyotype and t(9;11); and ‘high’: all other forms). The Supplementary Table 6 gives the list of 149 cDNAs selected by Bullinger et al. (2004) for building a binary prognostic factor, 147 of which have corresponding expression levels in Supplementary Table 3. Our analysis utilizes one single learning sample of n = 116 patients, 68 of whom died during the study period. The IPC weights are derived from a simple Kaplan–Meier estimate Gˆ of the censoring survivor function. For one patient a very late event was observed and we restrict the IPC weight for this patient to a value of five. Missing values in the expression matrix of all 6283 cDNAs and 116 patients are imputed using k = 10 nearest neighbor averaging (Troyanskaya et al., 2001), as implemented in the

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

Predictive modeling is the primary domain of ensemble methods, especially in situations where the number of covariates is large relative to the number of (uncensored) observations. A typical application is the construction of novel tumor classification schemes based on gene expression profiling data. One representative of such investigations is the recently published study of Bullinger et al. (2004) on AML patients. The main focus of this study was on the differentiation of previously unknown tumor subclasses by means of genetic information. Here, we try to construct ‘black box’ predictors for the survival time of AML patients incorporating both clinical and genetic information. Although the random forest or boosting estimate of the regression function f may be arbitrarily complex, some insight into the nature of the regression relationship is necessary in order to compare the fitted model with subject matter knowledge. In our second application, random forest and boosting are applied to data of a well-analyzed study on node-positive breast cancer. The estimated flexible regression functions are compared with previously published findings.

Survival ensembles

363

Table 1. Benchmark experiments for the GBSG-2 data: median performance for 100 bootstrap samples for the weighted mean (M), recursive partitioning (RP), a linear model (LM), random forest and L 2 -boosting for censored data with componentwise least squares p+ p+ p+ p+

=0 = 10 = 50 = 100

M 0.311 0.311 0.311 0.311

RP 0.311 0.311 0.311 0.311

LM 0.291 0.321 0.423 0.647

Random Forest 0.293 0.296 0.305 0.308

L 2 -boosting 0.289 0.299 0.303 0.310

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

Fig. 1. AML data: mean-difference plots (top) and scatterplots (bottom) of observed and predicted log-survival time of random forest and L 2 -boosting for censored data. The radius of the circles is proportional to the IPC weights. The dashed horizontal line in the lower panels is the weighted mean (with IPC weights) of the log-survival times, i.e. the prediction without any knowledge of the covariates.

R package pamr (Hastie et al., 2004). In total, 62 patients with IPC weights greater than zero had complete observations for the clinical variables and are used in the sequel. Random forest for censored data with 10 covariates randomly selected in each node of M = 250 trees and L 2 -boosting for censored data with componentwise linear models and AIC-based stopping criterion ( Mˆ = 350) were trained using both the eight clinical variables and the information covered by the 147 expression levels ( p = 155). The fit of both learners is depicted in Figure 1 and indicates a reasonable agreement between observed and predicted (log)-survival times for both algorithms. Both candidate models are compared with the naive overall mean prediction using a benchmark experiment following Hothorn et al. (2005b). From the learning sample L, 100 bootstrap samples are

364

T. H OTHORN ET AL .

drawn and the performance measures of all candidate models, i.e. the empirical risk defined in terms of the IPC weights, are evaluated on the same sample of out-of-bootstrap observations in an unreplicated complete block design. The candidate models have been fitted on the same bootstrap samples. The benchmark experiments are performed conditional on the IPC weights based on all observations, since we are interested in a comparison between the candidate models only. In order to investigate whether the molecular information of the expression levels helps to predict the survival time, we study in addition the performance of both algorithms when faced with a learning sample consisting of the clinical variables only (cRF and cL2B with p = 8). The joint and marginal distributions of the performance measures evaluated on the out-of-bootstrap observations are displayed in Figure 2, with median out-ofbootstrap errors of 2.451 (mean), 2.382 random forest, and 1.769 L 2 -boosting. In general, the performance distributions of the five candidate models show a global difference (asymptotic p-value < 0.001, Friedman test). All pairwise multiple comparisons based on Friedman rank sums (Wilcoxon–Nemenyi– McDonald–Thompson, see Hollander and Wolfe, 1999, Chapter 7.3) indicate that the naive prediction of the weighted mean is outperformed by AIC-based L 2 -boosting (adjusted p-value < 0.001). There is no evidence that the performance distributions of random forest and the weighted mean differ (adjusted p-value = 0.491). However, the distribution of the empirical risk of both ensemble methods is lower when only the eight clinical covariates are used (all adjusted p-values < 0.001). This supports the hypothesis that the raw gene expression levels do not improve the prediction of survival time. Bullinger et al. (2004) argue that the ‘likelihood and the duration of survival are likely to be fairly crude surrogates for the underlying biologic characteristics distinguishing prognostically relevant tumor subclasses’ and therefore propose an alternative strategy utilizing a prognostic variable obtained from a mix of cluster analysis and binary classification. 4.2

Node-positive breast cancer

A prospective, controlled clinical trial on the treatment of node-positive breast cancer patients was conducted by the German Breast Cancer Study Group (GBSG-2). A detailed description of the study is given in Schumacher et al. (1994). Patients not older than 65 years, with positive regional lymph nodes but no

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

Fig. 2. AML data: parallel coordinate plot and boxplots of the joint and marginal distributions of the L 2 risk with IPC weights evaluated on 100 out-of-bootstrap samples for the simple weighted mean (M), random forest and L 2 -boosting for censored data with componentwise least squares L 2 -boosting. In addition, the bootstrap errors for both ensemble methods based on the learning sample of the eight clinical covariates only are given (cRF and cL2B).

Survival ensembles

365

Fig. 3. GBSG-2 data: mean-difference plots of observed and predicted log-recurrence-free survival for all four candidate methods. The radius of the circles is proportional to the IPC weights.

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

distant metastases, were included in the study. Complete data on p = 7 prognostic factors for n = 686 women are used in Sauerbrei and Royston (1999) for prognostic modeling by means of multivariate fractional polynomials, i.e. flexible linear regression models based on transformed covariates. These findings will serve as the basis for the assessment of the diagnostic capabilities of survival ensembles. Observed hypothetical prognostic factors are age, menopausal status, tumor size, tumor grade, number of positive lymph nodes, progesterone receptor, estrogen receptor, and whether or not a hormonal therapy was administered. The recurrence-free survival time is the response variable of interest. The data are available in the R package ipred (Peters et al., 2002) and the IPC weights are derived from a simple Kaplan–Meier estimate Gˆ of the censoring survivor function. The weights are truncated to a maximal value of five for three very late events. The performance of four candidate algorithms is investigated: an ordinary linear model, fitted via IPC-weighted least squares; regression trees based on the IPC weights, as suggested by Molinaro et al. (2004) using the implementation in the R package rpart (Therneau and Atkinson, 1997); random forest for censored data (with five covariates randomly selected in each node of 100 trees); and L 2 -boosting for censored data, with componentwise linear models and AIC-based stopping criterion. The AIC-criterion for L 2 -boosting suggests stopping after the 86th boosting iteration. Figure 3 depicts a mean-difference plot of observed and predicted logarithms of recurrence-free survival for all four models. The figure leads to the impression that the relationship between the covariates and the recurrence-free survival time is relatively weak, a finding supported by an analysis with the Brier score in Hothorn et al. (2004).

366

T. H OTHORN ET AL .

Fig. 4. GBSG-2 data: joint and marginal distributions of the L 2 risk evaluated on 100 out-of-bootstrap samples for the weighted mean (M), random forest L 2 -boosting for censored data with componentwise least squares L 2 -boosting, recursive partitioning (RP), and a simple linear model (LM), for a number of additional random covariates p+ .

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

The performance of the four candidate models is compared by means of a benchmark experiment utilizing the framework given by Hothorn et al. (2005b) and described in Section 4.1. In order to study the stability of the models in high-dimensional situations, we chose a strategy in-between an analysis of the original data and a simulation experiment. We add p+ = (10, 50, 100) uncorrelated covariates, drawn from a uniform distribution, to the observed learning sample L and evaluate the performance using the out-of-bootstrap observations as described earlier. The results are depicted in Figure 4. Many-to-one comparisons based on Friedman rank sums indicate that, for the learning sample with only the original covariates ( p+ = 0), the linear model, boosting, and random forest perform better than the weighted mean (all adjusted p-values < 0.001). There is no evidence that the performance distributions of regression trees and the weighted mean differ (adjusted p-value = 0.965). Again, the relative improvement compared with the weighted mean is relatively small. For an increasing number of random covariates the linear model is heavily affected by overfitting but the ensemble methods are rather stable. For p+ = 50 additional random covariates, the bootstrap test set error of random forest and boosting is smaller than that for the weighted mean (both adjusted p-values = 0.001). However, there is only weak evidence that random forest for censored data performs better than the weighted mean for learning samples with p+ = 100 additional random covariates (adjusted p-value = 0.03); boosting cannot outperform the mean (adjusted p-value = 0.583) in this situation. The relative stability of regression trees is caused by the fact that the

Survival ensembles

367

Fig. 5. GBSG-2 data: scatterplots of selected covariates and predicted log-recurrence-free survival time obtained from random forest for censored data. A smoothing spline with 4 degrees of freedom is plotted. The radius of the circles is proportional to the IPC weights.

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

trees are pruned back to stumps or the root node most of the time. Figure 3 indeed shows (lower left panel) that only two predicted values are obtained from the tree stumps. Sauerbrei and Royston (1999) provide an in-depth analysis of the GBSG-2 data, focusing on fractional polynomials as interpretable but flexible regression models. They identify a nonlinear influence of the number of positive nodes, age, and progesterone receptor by visualization of the covariates and the corresponding (partial) linear predictions. We compare the estimated regression function fˆ provided by random forest and boosting with the findings reported in their paper by plotting the covariates against the predictions in Figures 5 and 6. These strategies were also applied for classification problems by Breiman (2001b) and Garczarek and Weihs (2003). The predicted log-recurrence-free survival time decreases with increasing number of positive lymph nodes (up to about 15 positive lymph nodes) for both random forest and boosting, in a nearly identical manner to the finding reported in Sauerbrei and Royston (1999). Both boosting and random forest suggest a relationship between age and survival time, namely, a decreasing risk for ages of 40–45 years and a nearly constant risk for older women, as in Sauerbrei and Royston (1999). A strong influence of the estrogen receptor is indicated by both ensemble methods. However, estrogen receptor measurements were not included in any of the models studied by Sauerbrei and Royston (1999). Progesterone receptor values (restricted to values less than 100 fmol/l) indicate a relationship with recurrence-free survival: very small values (less than about 10, say) are associated with short recurrence-free survival times, whereas

368

T. H OTHORN ET AL .

higher values indicate longer recurrence-free survival times. A similar finding is reported by Sauerbrei and Royston (1999). 5. D ISCUSSION The two algorithms presented in this paper extend ensemble prediction to censored data problems. Ensemble techniques have been developed in the past decade at the borderline between machine learning and statistics. Previous attempts to apply the main ideas to survival time data were bound to established key ingredients, such as the partial likelihood (Ridgeway, 1999), the Brier score for censored data (Benner, 2002), or survival trees (Hothorn et al., 2004; Ishwaran et al., 2004) and, consequently, inherited the associated difficulties. The general estimation framework of van der Laan and Robins (2003) and van der Laan and Dudoit (2003) allows for a sound theoretical formulation of the underlying risk optimization problems which can be solved with the new algorithms. Moreover, the framework enables us to apply well-known crossvalidation techniques for model evaluation (Keles¸ et al., 2004). Both ensemble algorithms are generic in the sense that arbitrary loss functions, e.g. absolute loss for any monotone transformation of the survival time T (including the identity), and other base learners can be implemented easily. For example, Henderson et al. (2001) investigate the prediction accuracy of survival models using a very intuitive loss function, having loss zero if the predicted survival time Tˆ satisfies 12 T  Tˆ  2T and loss one otherwise.

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

Fig. 6. GBSG-2 data: scatterpots of selected covariates and predicted log-recurrence-free survival time obtained from L 2 -boosting for censored data with componentwise least squares. A smoothing spline with 4 degrees of freedom is plotted. The radius of the circles is proportional to the IPC weights.

Survival ensembles

369

On the log-scale, this loss function reads L Henderson et al. (Y, ψ(X)) =

0, if |Y − ψ(X)|  δ, 1, if |Y − ψ(X)| > δ,

and the corresponding empirical risk can be minimized with respect to ψ by using pseudo-responses Y˜i − fˆm (Xi ), if |Y˜ − fˆm (X)|  δ, ˜ Ui = δ · sign(Y˜i − fˆm (Xi )), if |Y˜ − fˆm (X)| > δ, in Step 2 of the generic gradient boosting algorithm for censored data. For the GBSG-2 data, the meandifference plots and a comparison of observed versus predicted log-recurrence-free survival time for L 2 boosting optimizing the Huber and quadratic loss functions are shown in Figure 7. In fact, the fit based on the Huber loss function seems to be less variable compared to the fit obtained by optimizing quadratic loss. It should be noted that our implementations do not require an external choice of hyper-parameters, e.g. the number of boosting iterations. Another important property is that the random forest and the boosting algorithm reduce to their original full data form in the absence of censoring. In the uncensored data situation, the flexibility and stability of both the random forest and the boosting approach have been demonstrated in many benchmark experiments; we therefore restricted ourselves to a semiartificial benchmark experiment with varying number of covariates based on the GBSG-2 data. The main focus of our analysis of the AML and the GBSG-2 data is on the practical advantages of the methodology in terms of prediction accuracy and diagnostic ability. The results of flexible diagnostic modeling with fractional polynomials published by Sauerbrei and Royston (1999) could be reproduced for the GBSG-2 data. Thus, ensemble techniques are not just superb ‘black boxes’ in terms of prediction accuracy but can be used to investigate the nature of the regression relationship inherent in the data. We depicted simple partial relationships between one covariate and the predicted survival times. More advanced approaches for the visualization of complex regression relationships (Nason et al., 2004) are also applicable. When all observations are censored after a fixed time point, such as in a clinical trial running for a predefined period, the assumption P(C > T |X) > 0 stated in Section 2.2 is violated. Clearly, estimating the mean of a survival time is impossible when we never observe events after the end of follow-up. Here, we need to restrict our models to a truncated survival time as the response variable and to keep in mind that we can only derive conclusions about the part of the survival distribution for which we actually gathered information. The definition of the observed data loss function is the basis of all subsequent calculations. For the analysis of the AML and the GBSG-2 data, we used IPC weights obtained from a Kaplan–Meier estimate Gˆ of the censoring survivor function, i.e. an estimate based on T˜i and 1−i for observations i = 1, . . . , n (it should be noted that similar weighting schemes have been used to define measures of prediction accuracy for censored data, for example the Brier score by Graf et al., 1999). Molinaro et al. (2004) applied a Cox model to estimate the weights. This allows for modeling the censoring survivor function based

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

here with δ = log(2). Boosting based on this loss function is not possible because of nonconvexity— the gradient would be either zero or infinity. The squared error and also the Huber loss, suggested in the context of M-estimation and defined below, are convex functions with respect to the second argument. The Huber loss is a closer convex approximation (when properly scaled) of L Henderson et al. than the squared error loss. Specifically, 1 2 if |Y − ψ(X)|  δ, 2 (Y − ψ(X)) , L Huber (Y, ψ(X)) = δ(|Y − ψ(X)| − δ/2), if |Y − ψ(X)| > δ,

370

T. H OTHORN ET AL .

on information covered by a subset of the covariates. Robustness properties are studied theoretically in van der Laan and Robins (2003) and lead to DR-IPC weights as an alternative scheme. However, the practical implications of a misspecification of the weights, for example by omitting an important covariate when estimating the censoring distribution, and the advantages or disadvantages of parametric, semiparametric, or nonparametric modeling strategies need to be investigated by means of artificial simulation experiments. Another idea is to stabilize the estimate of the censoring distribution, and thus to stabilize the weights, by some form of ensemble technique prior to modeling or even simultaneously with the estimation of the regression function. Those issues are to be addressed in future research. 6. S OFTWARE All analyses were performed within the R system for statistical computing (R Development Core Team, 2004), version 2.0.1. A preliminary implementation of random forest for censored data is part of the R package party (Hothorn et al., 2005a). Until published on CRAN, implementations of the algorithms applied here are available from the authors upon request. ACKNOWLEDGMENTS We would like to thank two anonymous referees for their valuable comments. The work of T. Hothorn was supported by Deutsche Forschungsgemeinschaft under grant HO 3242/1-1, and S. Dudoit received support from the National Institutes of Health, grant NIH ROI LM07609.

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

Fig. 7. GBSG-2 data: mean-difference plots and scatterplots of observed and predicted log-recurrence-free survival for L 2 -boosting with Huber loss function and quadratic loss. The radius of the circles is proportional to the IPC weights.

Survival ensembles

371

R EFERENCES A LTMAN , D. G. AND ROYSTON , P. (2000). What do we mean by validating a prognostic model? Statistics in Medicine 19, 453–473. A MIT, Y. AND G EMAN , D. (1997). Shape quantization and recognition with randomized trees. Neural Computation 9, 1545–1588. BAIR , E. AND T IBSHIRANI , R. (2004). Semi-supervised methods to predict patient survival from gene expression data. PLoS Biology 2, 0511–0522. http://www.plosbiology.org. B EDRICK , E. J., E XUZIDES , A., J OHNSON , W. O. AND T HURMOND , M. C. (2002). Predictive influence in the accelerated failure time model. Biostatistics 3, 331–346.

B REIMAN , L. (1996). Bagging predictors. Machine Learning 24, 123–140. B REIMAN , L. (2001a). Random forests. Machine Learning 45, 5–32. B REIMAN , L. (2001b). Statistical modeling: the two cultures. Statistical Science 16, 199–231. B REIMAN , L. (2002). How to use Survival Forests. http://www.stat.berkeley.edu/users/breiman/. B RYAN , J., Y U , Z. AND VAN Biostatistics 5, 361–380. B UCKLEY, J.

AND JAMES ,

DER

L AAN , M. J. (2004). Analysis of longitudinal marginal structural models.

I. (1979). Linear regression with censored data. Biometrika 66, 429–436.

¨ , P. (2004). Bagging, boosting and ensemble methods. In Gentle, J. E., H¨ardle, W. and Mori, Y. (eds), B UHLMANN Handbook of Computational Statistics. Berlin: Springer-Verlag, pp. 877–907. ¨ B UHLMANN , P. (2006). Boosting for high-dimensional linear models. The Annals of Statistics 34 (in press). ¨ B UHLMANN , P. AND Y U , B. (2003). Boosting with L 2 loss: regression and classification. Journal of the American Statistical Association 98, 324–338. ¨ ¨ ¨ B ULLINGER , L., D OHNER , K., BAIR , E., F R OHLICH , S., S CHLENK , R. F., T IBSHIRANI , R., D OHNER , H. AND P OLLACK , J. R. (2004). Use of gene-expression profiling to identify prognostic subclasses in adult acute myloid leukemia. New England Journal of Medicine 350, 1605–1616. C OX , D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society, Series B 34, 187–202. D UDOIT, S. AND VAN DER L AAN , M. J. (2005). Asymptotics of cross-validated risk estimation in estimator selection and performance assessment. Statistical Methodology 2, 131–154. F RIEDMAN , J. H. (2001). Greedy function approximation: a gradient boosting machine. The Annals of Statistics 29, 1189–1202. G ARCZAREK , U. M. 18, 143–162.

AND

W EIHS , C. (2003). Standardizing the comparison of partitions. Computational Statistics

G RAF, E., S CHMOOR , C., S AUERBREI , W. AND S CHUMACHER , M. (1999). Assessment and comparison of prognostic classification schemes for survival data. Statistics in Medicine 18, 2529–2545. G RAY, R. J. (1992). Flexible methods for analyzing survival data using splines, with applications to breast cancer prognosis. Journal of the American Statistical Association 87, 942–951. H ASTIE , T. AND T IBSHIRANI , R. (2004). Efficient quadratic regularization for expression arrays. Biostatistics 5, 329–340. H ASTIE , T., T IBSHIRANI , R., NARASIMHAN , B. AND C HU , G. (2004). pamr: Prediction Analysis for Microarrays. R package version 1.24. http://CRAN.R-project.org.

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

B ENNER , A. (2002). Application of “aggregated classifiers” in survival time studies. In H¨ardle, W. and R¨onz, B. (eds), Proceedings in Computational Statistics: COMPSTAT 2002. Heidelberg: Physica-Verlag.

372

T. H OTHORN ET AL .

H ENDERSON , R. (1995). Problems and prediction in survival-data analysis. Statistics in Medicine 14, 161–184. H ENDERSON , R., J ONES , M. AND S TARE , J. (2001). Accuracy of point predictions in survival analysis. Statistics in Medicine 20, 3083–3096. H OLLANDER , M. AND W OLFE , D. A. (1999). Nonparametric Statistical Inference, 2nd edition. New York: John Wiley & Sons. H OTHORN , T., H ORNIK , K. AND Z EILEIS , A. (2005a). party: A Laboratory for Recursive Part(y)itioning. R package version 0.2-8. http://CRAN.R-project.org. ¨ , M. (2004). Bagging survival trees. Statistics H OTHORN , T., L AUSEN , B., B ENNER , A. AND R ADESPIEL -T R OGER in Medicine 23, 77–91.

H UANG , J. AND H ARRINGTON , D. (2005). Iterative partial least squares with right-censored data analysis: a comparison to other dimension reduction techniques. Biometrics 61, 17–24. I SHWARAN , H., B LACKSTONE , E. H., P OTHIER , C. E. AND L AUER , M. S. (2004). Relative risk forests for exercise heart rate recovery as a predictor of mortality. Journal of the American Statistical Association 99, 591–600. JAMES , I. (1998). Accelerated failure-time models. In Armitage, P. and Colton, T. (eds), Encyclopedia of Biostatistics. Chichester: John Wiley & Sons, pp. 26–30. K ELES¸ , S., VAN DER L AAN , M. J. AND D UDOIT, S. (2004). Asymptotically optimal model selection method for regression on censored outcomes. Bernoulli 10, 1011–1037. KOOPERBERG , C., S TONE , C. J. AND T RUONG , Y. K. (1996). Hazard regression. Journal of the American Statistical Association 90, 78–94. L E B LANC , M.

AND

C ROWLEY, J. (1999). Adaptive regression splines in the Cox model. Biometrics 55, 204–213.

L I , L. AND L I , H. (2004). Dimension reduction methods for microarrays with applications to censored survival data. Bioinformatics 20, 3406–3412. M OLINARO , A. M., D UDOIT, S. AND VAN DER L AAN , M. J. (2004). Tree-based multivariate regression and density estimation with right-censored data. Journal of Multivariate Analysis 90, 154–177. NASON , M., E MERSON , S. AND L E B LANC , M. (2004). CARTscans: a tool for visualizing complex models. Journal of Computational and Graphical Statistics 13, 807–825. O RBE , J., F ERREIRA , E.

AND

˜ -A NT ON ´ , V. (2003). Censored partial regression. Biostatistics 4, 109–121. N U´ NEZ

P ETERS , A., H OTHORN , T. AND L AUSEN , B. (2002). ipred: Improved predictors. R News 2, 33–36. http://CRAN. R-project.org/doc/Rnews/. R D EVELOPMENT C ORE T EAM . (2004). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org. R IDGEWAY, G. (1999). The state of boosting. Computing Science and Statistics 31, 172–181. R IPLEY, R. M., H ARRIS , A. L. AND TARASSENKO , L. (2004). Non-linear survival analysis using neural networks. Statistics in Medicine 23, 825–842. S AUERBREI , W. (1999). The use of resampling methods to simplify regression models in medical statistics. Journal of the Royal Statistical Society, Series C 48, 313–329. S AUERBREI , W. AND ROYSTON , P. (1999). Building multivariable prognostic and diagnostic models: transformation of the predictors by using fractional polynomials. Journal of the Royal Statistical Society, Series A 162, 71–94. S CHEMPER , M. (2003). Predictive accuracy and explained variation. Statistics in Medicine 22, 2299–2308.

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

H OTHORN , T., L EISCH , F., Z EILEIS , A. AND H ORNIK , K. (2005b). The design and analysis of benchmark experiments. Journal of Computational and Graphical Statistics 14, 675–699.

Survival ensembles

373

¨ S CHUMACHER , M., BASERT, G., B OJAR , H., H UBNER , K., O LSCHEWSKI , M., S AUERBREI , W., S CHMOOR , C., B EYERLE , C., N EUMANN , R. L. A. AND R AUSCHECKER , H. F., FOR THE G ERMAN B REAST C ANCER S TUDY G ROUP. (1994). Randomized 2 × 2 trial evaluating hormonal treatment and the duration of chemotherapy in node-positive breast cancer patients. Journal of Clinical Oncology 12, 2086–2093. S INISI , S. E. AND VAN DER L AAN , M. J. (2004). Deletion/substitution/addition algorithm in learning with applications in genomics. Statistical Applications in Genetics and Molecular Biology 3, Article 18. S TUTE , W. (1999). Nonlinear censored regression. Statistica Sinica 9, 1089–1102. T HERNEAU , T. M. AND ATKINSON , E. J. (1997). An introduction to recursive partitioning using the rpart routine. Technical Report 61, Section of Biostatistics, Mayo Clinic, Rochester. http://www.mayo.edu/hsr/techrpt/61.pdf. T HERNEAU , T. M. AND G RAMBSCH , P. M. (2000). Modeling Survival Data: Extending the Cox Model. New York: Springer.

L AAN , M. J. AND D UDOIT, S. (2003). Unified cross-validation methodology for selection among estimators: finite sample results, asymptotic optimality, and applications. Technical Report 130, Division of Biostatistics, University of California, Berkeley, California. http://www.bepress.com/ucbbiostat/paper130.

VAN DER

L AAN , M. J., D UDOIT, S. AND VAN DER VAART, A. W. (2004). The cross-validated adaptive epsilonnet estimator. Technical Report 142, Division of Biostatistics, University of California, Berkeley, California. http://www.bepress.com/ucbbiostat/paper142.

VAN DER

L AAN , M. J. AND ROBINS , J. M. (2003). Unified Methods for Censored Longitudinal Data and Causality. New York: Springer.

VAN DER

[Received April 4, 2005; first revision September 13, 2005; second revision November 4, 2005; accepted for publication December 7, 2005]

Downloaded from biostatistics.oxfordjournals.org at National Institute of Education Library, Serials Unit on April 20, 2011

T ROYANSKAYA , O., C ANTOR , M., S HERLOCK , G., B ROWN , P., H ASTIE , T., T IBSHIRANI , R., B OTSTEIN , D. AND A LTMAN , R. B. (2001). Missing value estimation methods for DNA microarrays. Bioinformatics 17, 520–525.