Iterative Bayesian Model Averaging: a method for the application of survival analysis to highdimensional microarray data
 Amalia Annest^{1},
 Roger E Bumgarner^{2},
 Adrian E Raftery^{3} and
 Ka Yee Yeung^{2}Email author
DOI: 10.1186/147121051072
© Annest et al; licensee BioMed Central Ltd. 2009
Received: 21 November 2008
Accepted: 26 February 2009
Published: 26 February 2009
Abstract
Background
Microarray technology is increasingly used to identify potential biomarkers for cancer prognostics and diagnostics. Previously, we have developed the iterative Bayesian Model Averaging (BMA) algorithm for use in classification. Here, we extend the iterative BMA algorithm for application to survival analysis on highdimensional microarray data. The main goal in applying survival analysis to microarray data is to determine a highly predictive model of patients' time to event (such as death, relapse, or metastasis) using a small number of selected genes. Our multivariate procedure combines the effectiveness of multiple contending models by calculating the weighted average of their posterior probability distributions. Our results demonstrate that our iterative BMA algorithm for survival analysis achieves high prediction accuracy while consistently selecting a small and costeffective number of predictor genes.
Results
We applied the iterative BMA algorithm to two cancer datasets: breast cancer and diffuse large Bcell lymphoma (DLBCL) data. On the breast cancer data, the algorithm selected a total of 15 predictor genes across 84 contending models from the training data. The maximum likelihood estimates of the selected genes and the posterior probabilities of the selected models from the training data were used to divide patients in the test (or validation) dataset into high and lowrisk categories. Using the genes and models determined from the training data, we assigned patients from the test data into highly distinct risk groups (as indicated by a pvalue of 7.26e05 from the logrank test). Moreover, we achieved comparable results using only the 5 top selected genes with 100% posterior probabilities. On the DLBCL data, our iterative BMA procedure selected a total of 25 genes across 3 contending models from the training data. Once again, we assigned the patients in the validation set to significantly distinct risk groups (pvalue = 0.00139).
Conclusion
The strength of the iterative BMA algorithm for survival analysis lies in its ability to account for model uncertainty. The results from this study demonstrate that our procedure selects a small number of genes while eclipsing other methods in predictive performance, making it a highly accurate and costeffective prognostic tool in the clinical setting.
Background
Introduction and Previous Work
Until recently, oncologists relied primarily on tumor stage and morphology to help outline an appropriate course of treatment for their cancer patients. Malignant tumors were generally resected in operable cases, and followup radiation therapy was provided to victims exhibiting advancedstage diseases. This methodology proved problematic in that a number of lowrisk patients experienced cancer recurrence or death within a short time frame, while a contingent of highrisk patients went into permanent remission despite the bleak nature of their original prognoses. This indicated a need to explore other indicators by which doctors could understand the underlying prognosis of a given disease and decide on a treatment plan that would optimize the patient's chances for survival.
Microarray technology provides a promising avenue. The availability of thousands of gene expression levels has enabled the pursuit of a new direction in cancer research. In particular, gene expression patterns can be thought of as multidimensional quantitative "expression phenotypes" which can in turn be correlated with clinical outcome. Because a single microarray can measure the expression levels of tens of thousands of genes simultaneously, the challenge lies in the development of data mining methods and tools to extract biological meaning from this immense amount of data. More specifically, the aim is to filter the expression dataset down to the smallest possible subset of accurate predictor genes. Reducing the number of predictor genes both decreases clinical costs and mitigates the possibility of overfitting due to high intervariable correlations [1].
The most common approach to identify a manageable group of predictor genes is called feature selection, in which a subset of relevant "features" (or variables) is selected from the full dataset in order to produce a robust learning model [2, 3]. A welldesigned feature selection algorithm will choose a small set of variables that is highly predictive of clinical outcome. Univariate feature selection methods evaluate the usefulness of each variable on an individual basis. Examples of univariate techniques include the ttest [4], the signaltonoise ratio [5], the Cox proportional hazards model [6], threshold number of misclassification (TNoM) score [7], the betweengroups to withingroups sum of squares (BSS/WSS) ratio [8], and mean aggregate relevance [9]. Multivariate methods are more sophisticated in that they perform combinatorial searches within the feature subspace to evaluate the effectiveness of groups of genes. Examples include Recursive Feature Elimination (RFE) [10], genetic algorithms [11–13], floating search [14], and topscoring pair methods [15, 16]. Despite some evidence to the contrary (e.g., [17]), multivariate selection algorithms are generally preferable to univariate ones because they cut down on dependencies between variables and often lead to models with fewer predictive variables [18, 19]. However, selecting multivariate features from microarray data is nontrivial since the number of patient samples is often limited (usually under a hundred) and the number of genes is large (usually tens of thousands).
Subsequent to or concurrent with the feature selection process, a supervised machine learning technique can be applied to generate a predictive function using the selected variables from a set of training data [20, 21]. In a supervised learning algorithm, the input is a set of training samples paired with the corresponding labels of those samples. If the labels are exhaustive discrete classes to which the samples belong (e.g. "survived beyond five years" and "died before five years"), then the learning model is a classifier (for a review of classification techniques in supervised machine learning, see [22]). With microarray data, the most common approach is to apply a classification algorithm in which the patients are split into subcategories corresponding to different prognoses or diagnoses. In general, the subcategories are static and based on thresholds associated with some clinical variable (e.g., time to metastases). Classification studies on microarray data have used gene expression levels to distinguish diseased tissue samples from normal ones [23, 24], identify cancer subtypes [5, 25], and assign discrete risk groups for survival prognosis [13, 19, 26–28]. See Hu et al. [29] for a comparative analysis of classification methods for microarray data. Depending on whether the classifier is used to select relevant features, feature selection methods can also be divided into filter and wrapper methods [30]. Wrapper methods utilize the classifiers as evaluation functions and search for the optimal gene set for classification. In contrast, filter methods rely on general characteristics of the training data to select genes without involving any classifier for evaluation. Many filter methods evaluate a gene based on its discriminative power for the target classes without considering its correlations with other genes. Wrapper algorithms can perform better than filter algorithms, but they typically require orders of magnitude more computation time.
In cancer research, gene expression data is often reported in tandem with time to event information (such as time to metastasis, death, or relapse). In order to take advantage of these continuous clinical variables under a supervised framework, survival analysis can be applied. Survival analysis on microarray data differs from classification in that the sample labels are continuous rather than discrete. The overall goal in survival analysis research is to create the strongest predictive model of patient survival, and the most important components of this process are feature selection and model construction. In the context of survival analysis, a model refers to a set of selected genes whose regression coefficients have been calculated for use in predicting survival prognosis [31]. In the application of survival analysis to highdimensional microarray data, a feature selection algorithm identifies this subset of genes from the gene expression training dataset. These genes are then used to build a statistical model for the continuous time to event data [32]. The choice of feature selection algorithm determines which genes are chosen and the number of predictor genes deemed to be relevant, whereas the statistical model gives the distribution of the time to the event.
In recent years, a number of studies have applied survival analysis to microarray data. Beer et al. [33] used univariate Cox proportional hazards regression along with leaveoneout cross validation on an 86sample lung cancer dataset to develop a risk index based on 50 genes that successfully divided an independent test set of patients into high and lowrisk groups. Lu et al. [34] improved on these results by using multivariate Cox proportional hazards model with bootstrap resampling and forward selection to obtain a 64gene model that yielded a greater predictive accuracy than Beer et al. A popular approach to deal with high dimensionality in survival analysis is dimension reduction. For example, Bair and Tibshirani [35] proposed a semisupervised version of principal components analysis that is capable of generating a continuous predictor of patient survival. Their algorithm consistently selected fewer than 20 genes and successfully divided patients into high and lowrisk groups in four different cancer subtypes: lymphoma, breast cancer, lung cancer, and acute myeloid leukemia. Partial least squares (PLS) reduces the dimension of the original variables by constructing a smaller collection of latent variables that are linear combinations of the original variables. The application of PLS in conjunction with the Cox proportional hazards model in survival analysis to microarray data has been investigated [36, 37]. A drawback of dimension reduction techniques is that usually a relatively large number of genes (variables) are selected in the reduced dimension space. Penalized methods such as LASSO (least absolute shrinkage and selection operator) [38, 39] is a variable selection method, and hence is an alternative to dimension reduction techniques as it can be used when the number of samples is smaller than the number of variables (genes). Zhang et al. [40] studied the theoretical properties of an adaptive LASSO method for the Cox proportional hazards model. Kaderali et al. [41] proposed a multivariate Cox regression model embedded in a Bayesian framework that combines dimension reduction and regression in one single step. They used a hierarchical prior distribution that is strongly peaked around zero on the regression parameters so as to produce a small number of relevant genes with nonzero regression parameters. A distinctive characteristic in Kaderali et al. is that they assume that the constant baseline hazard rate in the Cox proportional hazards model is known, and aim to directly predict survival times of patients. Recently, Bovelstad et al. [42] compared the prediction performance of seven methods that are based on the Cox proportional hazards model over three microarray datasets, and showed that ridge regression has the overall best performance. In their empirical studies, Bovelstad et al. focused on prediction accuracy instead of the number of selected genes.
The accelerated failure time (AFT) model [43] is a linear regression model in which the response variable is the logarithm or a known monotone transformation of event times. Unlike the Cox proportional hazards model which assumes that the ratio of the hazard functions does not depend on time and the baseline hazard is unspecified, one can directly predict event times using the AFT model and hence, is a useful alternative to the Cox model. However, the AFT model has not been widely used in practice due to difficulties in computing the regression parameters even when the number of variables is small. Recently, Huang and colleagues [44, 45] studied the use of penalized methods in the AFT model for survival analysis. However, these proposed methods were not applied to microarray datasets in which thousands or tens of thousands of variables are available. Datta et al. [46] investigated the performances of LASSO and PLS on microarray data using the AFT model. They showed that LASSO performed better than PLS when there are many noise variables in their simulation studies.
Our Contributions
A problem with most feature selection algorithms used to produce continuous predictors of patient survival is that they fail to account for model uncertainty. With thousands of genes and only tens to hundreds of samples, it often happens that a number of different models describe the data about equally well. In this paper, we apply the Bayesian Model Averaging (BMA) method [47, 48] to select a subset of genes for survival analysis on microarray data. Instead of choosing a single model and proceeding as if the data were actually generated from it, BMA combines the effectiveness of multiple models by taking the weighted average of their posterior distributions. In addition, BMA consistently identifies a small number of predictive genes [19, 31], and the posterior probabilities of the selected genes and models are available to facilitate an easily interpretable summary of the output. Yeung et al. [19] extended the BMA algorithm to classify highdimensional microarray; they dealt with the very large number of potential predictors using an iterative approach. Here we further extend their iterative BMA method to survival analysis. In particular, we developed and implemented the iterative BMA method for survival analysis as a Bioconductor package, and we also demonstrated our algorithm on two cancer datasets. Our results reveal that iterative BMA consistently selects a small number of predictor genes while providing greater predictive accuracy than other algorithms, and the models themselves are simple and amenable to biological interpretation.
Methods
Data
Breast Cancer
Summary of Breast Cancer and DLBCL Datasets
Dataset  Total Number of Samples  # Training Samples  # Validation Samples  Number of Genes 

Breast Cancer  295  61  234  4919 
DLBCL  240  160  80  7399 
Lymphoma
Our second dataset consists of tumor samples from 240 patients diagnosed with diffuse large BCell lymphoma (DLBCL) [50]. Roughly 60% of DLBCL victims who are treated with chemotherapy do not survive, and the disease comprises 30–40% of all nonHodgkin lymphomas [51, 52]. This DLBCL dataset was generated and first analyzed by Rosenwald et al. [50], and the expression profiles from 7399 genes along with corresponding patient information can be downloaded from their supplemental website http://llmpp.nih.gov/DLBCL/. The raw data were processed with "lymphochip" cDNA microarrays [53], which are specialized to include genes that are known to be preferentially expressed within the germinal centers of lymphoid organs. Survival times ranged from 0 to 21.8 years, with a median of 2.8 years across all samples. Of the 240 patients, only 102 (42.5%) were still alive at the final followup visit. Rosenwald et al. randomly divided the dataset into 160 training samples and 80 validation samples, and we have chosen to preserve their division in order to allow a direct comparison of results. See Table 1 for a summary of the breast cancer and DLBCL datasets.
Bayesian Model Averaging (BMA)
The strength of BMA lies in its ability to account for model uncertainty, an aspect of analysis that is largely ignored by traditional stepwise selection procedures [47]. These traditional methods tend to overestimate the goodnessoffit between model and data, and the model is subsequently unable to retain its predictive power when applied to independent datasets [31, 54]. BMA attempts to solve this problem by selecting a subset of all possible models and making statistical inferences using a weighted average of these models' posterior distributions.
There are three issues to consider before Equation (1) can be applied: obtaining the subset S of models to be included, estimating the value of Pr(Ψ  TD, M_{i}), and estimating the value of Pr(M_{i}  TD). The remainder of this subsection will address these issues.
One challenge with BMA is the sheer number of models that could potentially be explored by the algorithm, especially when dealing with microarray data. If there are G candidate explanatory genes in the expression set, then there are 2^{G} possible models to consider. When working with tens of thousands of genes, such an undertaking is computationally intractable. In order to discard the noncontributory models and obtain a subset that approximates an average over all 2^{G} possibilities, Raftery [47] proposed to use the regression by leaps and bounds algorithm from Furnival and Wilson [55]. This algorithm takes a userspecified input "nbest" and efficiently returns the top nbest models of each size (maximum 30 variables). Following application of the leaps and bounds algorithm, the Occam's window method of Madigan and Raftery [56] can be used to reduce the set of models. After identifying the strongest model returned by the leaps and bounds algorithm, the procedure can eliminate any model whose posterior probability is below the cutoff point in relation to the best model. The cutoff point can be varied, but the default is 20; that is, a model must be at least 1/20 as likely as the strongest model in order to be retained. Once this step is complete, the remaining group of models constitutes the set S to be used in Equation (1).
An exact calculation of the predictive distribution Pr(Ψ  TD, M_{i}) requires an integration over the vector of regression parameters θ_{i}:
Pr(Ψ  TD, M_{i}) = ∫ Pr(Ψ  θ_{i}, TD, M_{i}) Pr(θ_{i}  TD, M_{i}) dθ_{i}. (2)
While certain techniques such as the Markov Chain Monte Carlo (MCMC) methods have been used in survival analysis to obtain a more exact predictive distribution [57], the MLE requires fewer computational resources and has been deemed sufficient for the purpose of averaging over contending models [31, 58–60].
Finally, a calculation of the posterior probability of model M_{i} given the training data TD involves an integral whose value is impossible to evaluate exactly. Bayes' theorem yields Equation (4), which represents the posterior probability of model M_{i} given TD:
Pr(M_{i}  TD) ∝ Pr(TD  M_{i}) Pr(M_{i}), (4)
where
Pr(TD  M_{i}) = ∫ Pr(TD  θ_{i}, M_{i}) Pr(θ_{i}  M_{i}) dθ_{i}. (5)
In equation (6), n represents the number of records in the data, k_{i} is the number of regression parameters in model M_{i}, and O(1) is the error term. The approximation is more accurate for many practical purposes than its O(1) error term suggests for certain reasonable choices of the prior distribution Pr(θ_{i}  M_{i}) [47, 61]. Raftery [62] gave further empirical evidence for the accuracy of this approximation. This method is implemented in the bic.surv function that is part of the BMA R package available at http://cran.rproject.org/web/packages/BMA/index.html
The posterior probability of gene x_{i} is the sum of the posterior probabilities of all models in the subset S that include gene x_{i}.
BMA for Survival Analysis
Volinsky et al. [31] applied the Bayesian Model Averaging methods [47, 48] to survival analysis. They assessed a patient's risk of stroke by using BMA to select variables in Cox Proportional Hazards models [6]. The data were made available by the Cardiovascular Health Study and included 23 variables (e.g., age, smoking history, and blood pressure) that may contribute to a patient's chances of experiencing a stroke. BMA selected a total of 5 models and 11 predictive variables, including diuretic, aspirin use, diabetes, stenosis, and timed walk. Patient risk scores were calculated by taking the weighted average of the risk scores for each of the top five contending models. The patients were then assigned to either the highrisk, mediumrisk, or lowrisk group based on the empirical 33^{rd} and 66^{th} percentile cutoff points in the risk scores of the training set. To assess performance, Volinsky et al. [31] created an analogue to the logscore called the partial predictive score (PPS). The PPS for BMA was compared against the PPS for the top BMA model (that is, the single model of the top five BMA models with the highest posterior probability) and against the PPS of the model returned by stepwise backward elimination. BMA exhibited the highest PPS, with a prediction mechanism 15% more effective than the top model alone and 3.5% more effective than the stepwise procedure. Furthermore, the patients assigned to a risk group using BMA experienced fewer strokes in the lowrisk group and more strokes in the highrisk group when compared with the other two methods.
Extending BMA for HighDimensional Microarray Data
Iterative BMA for Classification
The BMA implementation described above is incompatible with microarray data. This is because the typical microarray dataset contains thousands or even tens of thousands of genes, but the leaps and bounds algorithm from Furnival and Wilson [55] tends to become slow when there are more than 45 variables or so. One common solution is to use stepwise backward elimination to reduce the number of genes down to 30, but this is not applicable in a situation where the number of predictive variables is greater than the number of samples. Yeung et al. [19] developed an iterative BMA algorithm that takes a rankordered list of genes and successively applies the traditional BMA algorithm until all genes up to a userspecified value p (G_{1}, G_{2}, ..., G_{p}) have been processed. The authors begin by using the ratio of betweengroup to withingroup sum of squares (BSS/WSS) [8] to rankorder the genes from the microarray dataset. As the algorithm iterates, genes with a high posterior probability (equation (7)) are retained while genes with a low posterior probability are eliminated. The default threshold for inclusion is set to 1%; genes whose posterior probabilities are less than 1% are discarded.
Iterative BMA for Survival Analysis
In this article, we report our efforts in extending the iterative BMA method to survival analysis, which include a number of modifications to the algorithm. First, instead of applying the BSS/WSS technique [8] to rankorder the genes in the preprocessing step, we use the Cox Proportional Hazards Model [6] to rank each individual gene. Cox regression is a popular choice in the realm of survival analysis due to its broad applicability and capacity for handling censored data. It is a semiparametric method that quantifies the hazard rate for a subject s at time T as follows:
λ(T  p_{s}) = λ_{0}(T)exp(p_{s}θ). (8)
In equation (9), R_{s} is the risk set at time t_{s} (where the risk set consists of individuals who have not yet experienced the event of interest), and δ_{i} is an indicator for whether subject i is censored. Once the regression parameters in the Cox model are estimated by maximizing the partial likelihood, the genes can be ranked in descending order of their log likelihood.
Furthermore, we have incorporated an additional heuristic in which the models that are discarded due to the adaptive threshold are reconsidered when all the iterations are completed. Specifically, we applied the Occam's window method of Madigan and Raftery [56] to reduce the set of models remaining from the last iteration of bic.surv and the models that are discarded due to the adaptive threshold. We also recomputed the posterior probabilities of the models and the genes accordingly. In our software implementation, this heuristic is available as an option called "keepRmModels" and this option is set to FALSE by default.
Assessment
In equation (10), ${\widehat{\theta}}_{i}$ represents the vector of regression parameters for model M_{ i }and ${x}_{j}^{v}$ refers to the expression score of each gene x_{ j }within model M_{ i }for a patient in the validation dataset. Therefore, the risk score is computed by multiplying the expression scores of all genes included in model M_{ i }by their corresponding predictor coefficients, adding these x_{ j }b_{ j }terms together, weighing this number by the posterior probability of each model M_{ i }and summing over all contending models in the set S. Note that the predictor coefficients and the model posterior probabilities are all determined from the training data. Our implementation employs uses the userspecified "cutPoint" for defining high versus lowrisk groups (e.g., a cutPoint of 60 means the lower 60% of scores will be deemed lowrisk, and the upper 40% will comprise the highrisk group) using the risk scores of patients in the training data.
The KaplanMeier survival curves [63], in which the proportions of surviving patients in each risk group are plotted against successive time intervals, are used to illustrate our results. An advantage of the KaplanMeier curve is that it takes censored data into consideration: small vertical tickmarks represent losses where patient data were censored. In addition, we measured predictive performance with the pvalue calculated from the logrank test using the central chisquare distribution. The logrank test calculates a pvalue testing the null hypothesis that the survival curves from the high and lowrisk groups are identical. Therefore, a significant pvalue indicates that the two risk groups are distinct.
Selection of Input Parameters
10run/10fold cross validation results on the DLBCL dataset for cutPoint = 60 and maxNvar = 25.
p (# genes)  nbest  Average pvalue  pvalue stdev  Average chisquare value  chisquare stdev 

500  10  0.385  0.308  2.048  2.831 
500  20  0.398  0.313  1.853  2.349 
500  50  0.329  0.291  2.211  2.842 
500  100  0.320  0.294  2.414  2.735 
1000  10  0.313  0.303  2.648  3.139 
1000  20  0.369  0.308  2.107  2.588 
1000  50  0.307  0.303  2.958  3.251 
1000  100  0.310  0.271  2.493  3.040 
For the results shown in the rest of this paper, we adopted the optimal input parameters (p = 1000 and nbest = 50) determined from the 10run/10fold cross validation study on the DLBCL dataset shown in Table 2 in addition to the chosen 60% cutPoint to define risk groups. On the DLBCL data, we used the default BMA window size (represented by the input parameter maxNvar) of 25. However, on the breast cancer data, we were unable to use the default maxNvar value of 25. During preliminary analyses (data not shown), we found that training sets with fewer than 100 samples tended to result in fatal errors at higher values of maxNvar. These errors occur because smaller matrices lead to instability in fitting the data. These singularity errors can be largely mitigated by reducing the number of variables in the active BMA window. Since the breast cancer training set is relatively small (61 samples vs. 160 samples in the DLBCL data), we reduced the value of maxNvar from 25 to 15 variables in order to avoid convergence errors caused by matrix singularity.
Results and discussion
Breast Cancer Data
Genes selected by the iterative BMA algorithm and their corresponding posterior probabilities, univariate log likelihood rankings, and descriptions on the breast cancer data (p = 1000, nbest = 50, maxNvar = 15, and cutPoint = 60).
Selected genes  Posterior Probability (%)  Univariate Cox ranking  Gene description 

NM_000767  100.0  437  cytochrome P450, subfamily IIB (phenobarbitalinducible) 
NM_002019  100.0  533  fmsrelated tyrosine kinase 1 (vascular endothelial growth factor/vascular permeability factor receptor) 
Contig47102_RC  100.0  564  no description available 
NM_013989  100.0  765  deiodinase, iodothyronine, type II (DIO2), transcript variant 1, mRNA 
NM_018965  100.0  935  triggering receptor expressed on myeloid cells 2 
NM_021151  99.0  956  carnitine O octanoyltransferase 
AF063936  43.9  984  putative neuronal cell adhesion molecule 
NM_004911  40.6  998  protein disulfide isomerase related protein (calciumbinding protein, intestinalrelated) 
NM_014862  29.5  994  KIAA0307 gene product 
Contig40146  19.0  996  wi84e 12.x1 NCI_CGAP_Kid12 Homo sapiens cDNA clone IMAGE : 2400046 3' similar to SW: RASD_DICDI P03967 RAS LIKE PROTEIN RASD;, mRNA sequence 
NM_012319  17.6  993  LIV1 protein, estrogen regulated 
NM_002411  13.1  995  secretoglobin, family 2A, member 2 (SCGB2A2), mRNA 
NM_003645  10.7  997  fattyacidCoenzyme A ligase, very longchain 1 
NM_012415  10.3  1000  RAD54 homolog B (S. cerevisiae), transcript variant 1, mRNA (cDNA Clone, ORF Clone) 
NM_015972  9.4  999  polymerase (RNA) I polypeptide D, 16 kDa 
The number of censored and uncensored breast cancer patient samples in each risk group, along with the total number of censored and uncensored patients and the total number of patients in the high and lowrisk categories.
Censored  Uncensored  Total  

High risk  69  38  107 
Low risk  110  17  127 
Total  179  55 
Furthermore, we investigated the effect of the adaptive threshold heuristic in which the posterior probability threshold is increased temporarily to ensure that at least one gene is removed in each iteration. Our exploration showed that the adaptive threshold played an important role in the results shown in Table 3. Specifically, many of the high univariate ranked genes were removed due to the adaptive threshold. If we turned off the adaptive threshold heuristic, there will only be a single iteration of bic.surv since all the top univariate ranked genes produced posterior probabilities greater than the 1% threshold. Without the adaptive threshold, our analysis showed that using the top 15 univariate genes produced less distinct risk groups (pvalue = 6.21e4, chisquare = 11.712 for n = 234). We have also explored the keepRmModels=TRUE heuristic in which all the models discarded due to the adaptive threshold were reconsidered by the Occam's window method. In this case, our analysis yielded mostly highranked univariate genes (32 genes spanning across 217 models) that produced slightly less distinct risk groups (pvalue = 2.31e3, chisquare = 9.283 for n = 234). Please refer to Additional file 1 for more detailed results using this heuristic.
A comparison of the pvalues and chisquare statistics from the logrank test and numbers of genes selected across different survival analysis methods on the full breast cancer validation set of van de Vijver et al. [49] (n = 295), and the partial breast cancer validation set used in this work with 61 overlapping samples removed (n = 234).
heuristic  # genes  n = 234  n = 234  n = 295  n = 295  

pvalue  chisquare  pvalue  chisquare  
iterative BMA  adaptive threshold, keepRmModels = FALSE  15  7.264E05  15.714  3.382E10  39.441 
iterative BMA  top 5 genes with 100% posterior probabilities  5  9.063E06  19.699  1.143E10  41.559 
iterative BMA  adaptive threshold, keepRmModels = TRUE  32  2.312E03  9.283  9.875E08  28.398 
Bair & Tibshirani (2002) Principle Components  NA  5  3.280E03  8.645*  3.120E05  17.343* 
Method of van't Veer et al. (2002) (as calculated by Bair & Tibshirani)  NA  70  1.050E02  6.548*  9.340E03  6.757* 
DLBCL Data
Genes selected by the iterative BMA algorithm and their corresponding posterior probabilities, univariate log likelihood rankings, and descriptions on the DLBCL dataset (p = 1000, nbest = 50, maxNvar = 25, and cutPoint = 60).
Selected genes  Posterior Probability (%)  Univariate Cox ranking  Gene description 

BC012161  100.0  1  septin 1 
D42043  100.0  4  KIAA0084 protein 
X53505  100.0  41  ribosomal protein S12 
BF129543  100.0  49  ESTs, weakly similar to A47224 thyroxinebinding globulin precursor 
D13666  100.0  73  osteoblast specific factor 2 (fasciclin Ilike) 
M83664  100.0  93  MHC, class II, DP beta I 
AK000978  100.0  101  hypothetical protein FLJ10116 
AF009615  100.0  116  a disintegrin and metallo proteinase domain 10 
AK027711  100.0  123  hypothetical protein MGC3234 
LC_24015  100.0  129  no description available 
K01144  100.0  140  CD74 antigen (invariant polypeptide of MHC, class II antigen associated) 
U68418  100.0  181  branched chain aminotransferase 2 
D88532  100.0  213  phosphoinositide3kinase, regulatory subunit, polypeptide 3 (p55, gamma) 
NM_022551  100.0  223  ribosomal protein S18 
U18259  100.0  242  MHC, class II transactivator 
X64707  100.0  243  ribosomal protein L13 
NM_006312  100.0  278  nuclear receptor co repressor 2 
M58297  100.0  385  zinc finger protein 42 (myeloid specific retinoic acidresponsive) 
LC_26524  100.0  473  no description available 
AK022743  100.0  499  hypothetical protein FLJ12681 
NM_005347  100.0  518  heat shock 70 kDa protein 5 (glucoseregulated protein 78 kDa) 
D83492  100.0  632  EphB 
AK025754  100.0  652  HP1BP74 
AA747694  81.7  885  ESTs, weakly similar to ALU SUBF J 
U70981  9.2  1000  interleukin 13 receptor, alpha 2 
The number of censored and dead DLBCL patients in each risk group, along with the total number of censored and dead patients and the total number of patients in the high and lowrisk categories.
Censored  Died  Total  

High risk  3  21  24 
Low risk  27  29  56 
Total  30  50 
A comparison among three studies of the number of genes selected, the corresponding pvalues and chisquare statistics in survival analysis on the DLBCL dataset.
Number of Genes  pvalue  chisquare  

iterativeBMA  25  0.00139  10.221 
Bair & Tibshirani (2004)  17  0.00124  10.430* 
Rosenwald et al. (2002)  37–1333  0.009 – 0.11  2.554 – 6.823* 
As with the breast cancer data, we restricted the number of predictor variables in the DLBCL training set to include only those genes with posterior probabilities of 100%. Because 23 of the 25 genes originally selected by iterativeBMAsurv belonged in this category, the outcome was almost identical to the results presented in Figure 6. This information suggests that rerunning the algorithm with a reduced set of predictor variables may be more worthwhile when relatively fewer genes are calculated to have the maximum posterior probability.
Comparison with Other Methods
Here we compared our iterativeBMAsurv results to that from using ridge regression which was shown to produce the best overall prediction accuracy in a recent empirical study from Bovelstad et al. [42]. Specifically, Bovelstad et al. compared the prediction performance of seven methods that are based on the Cox proportional hazards model, including univariate selection, forward selection, principal components regression, supervised principal components regression, partial least squares regression, ridge regression and LASSO. Since Bovelstad et al. reported results produced by randomly splitting three microarray datasets into training and test sets, our results are not directly comparable with theirs. Therefore, we downloaded their software implementation written in Matlab from their supplementary web site [64], ran their software on the breast cancer and DLBCL datasets from which we derived our results, and then assessed the results using the logrank test. On the breast cancer data, ridge regression produced a pvalue of 0.00340 using all 4919 genes. In contrast, our iterativeBMAsurv algorithm produced a much more significant pvalue of 9.063e06 using only 5 genes. On the DLBCL data, ridge regression produced a pvalue of 0.000380 using all 7399 genes. In contrast, our iterativeBMAsurv algorithm produced a pvalue of 0.001389 using 25 genes. Bovelstad et al. focused on prediction accuracy instead of the number of selected genes. In summary, ridge regression may produce good prediction accuracy, but it is not a variable selection algorithm. On the other hand, our iterativeBMAsurv algorithm typically selects a small number of genes and produces good prediction accuracy.
Conclusion
A summary of the results from the application of the iterative BMA algorithm to the DLBCL dataset, the partial nonoverlapping breast cancer dataset (n = 234), and the full overlapping breast cancer dataset (n = 295).
Number of Genes  Number of Models  pvalue  chisquare  

DLBCL  25  3  1.389e03  10.221 
Breast Cancer n = 234 (15 genes)  15  84  7.264e05  15.714 
Breast Cancer n = 295 (15 genes)  15  84  3.382e10  39.441 
Breast Cancer n = 234 (5 genes)  5  2  9.063e06  19.699 
Breast Cancer n = 295 (5 genes)  5  2  1.143e10  41.559 
Our results showed that genes with poor univariate rankings are often selected by the iterative BMA algorithm using the adaptive threshold heuristic (see Tables 3 and 6). Our analysis showed that genes with dramatic difference in the univariate rankings may in fact have similar goodnessoffit (i.e. comparable log likelihood) when fitted to the univariate Cox proportional hazards model. On the breast cancer data, we showed that our selected set of genes (with poor univariate rankings) resulted in patients being assigned to more distinct risk groups than the top univariate genes. While it is true that the genes and models selected by the iterative BMA procedure are contingent upon the initial univariate rankings, all p topranked univariate genes may be included in the models selected by our iterativeBMAsurv method. Our results showed that setting the parameter p to a large value (e.g., 1000) generally yields high prediction accuracy. In addition to the parameter p, our algorithm requires the input of a few other userspecified parameters. One example is nbest, which is used by the leaps and bounds algorithm from Furnival and Wilson [55] to isolate the nbest strongest models. Higher values of nbest increase the computation time, but overly restrictive values undermine predictive accuracy by failing to return potentially contributory models. We found that nbest = 20, 50, and 100 generally yielded good results, with a value of 50 exhibiting the ideal tradeoff between predictive power and computational efficiency on the DLBCL data. For example, it takes about 1.5 hours to run the iterative BMA algorithm with p = 1000 genes and nbest = 50 on a machine with 2 gigabytes of RAM and a 2.0 GHz Intel dual core processor. Reducing nbest to 20 cuts the running time down to 40 minutes, but the pvalue and chisquare statistic representing the difference between risk groups is slightly less favorable. On the other hand, setting nbest to 100 significantly increases the computation time with no appreciable improvement in prediction. Cross validation can be used to determine the optimal input parameters for each dataset.
While the results obtained in this study are encouraging, the iterative BMA algorithm for survival analysis presents some limitations and areas for future development. The mathematical calculations conducted by the BMA methods are close approximations, but they could be computed with greater precision. For example, the maximum likelihood estimate of equation (3) provides a sufficient approximation to the predictive distribution, but the more computationally intensive Markov Chain Monte Carlo methods might yield the true predictive distribution with greater accuracy [31]. The approximation of the posterior model probabilities calculated in equation (6) could also be improved [62]. Another area for future work lies in the assessment of the performance from different computational methods. Currently, computational methods are evaluated by comparing the separation between different risk groups using the logrank test (e.g. [35, 50]). In our work, we divided patients into the high and low risk groups using cutPoint = 60% and evaluated computational methods based on the pvalues and chisquare statistics computed using the logrank test. Note that we are only comparing the pvalues over identical test sets. Comparing pvalues across different test sets could potentially be misleading, and hence, we recommend using both the pvalues and chisquare statistics when evaluating different computational methods. Further work is required to investigate the optimal number of risk groups and to propose statistical methods for the assessment of different computational methods across different test sets.
The iterative BMA algorithm for survival analysis is easy to use, computationally efficient, and highly accurate. It identifies a handful of predictor variables from vast amounts of microarray data, making it a costeffective diagnostic tool in the clinical setting. In terms of future work, we would like to collaborate with cancer biologists to validate the predictor genes selected by applying the iterativeBMAsurv algorithm to microarray data, and to assess the prediction accuracy of our methodology on PCR data generated using independent patient samples. Furthermore, we would like to extend the iterative BMA algorithm to other types of highthroughput data such as proteomics data produced from mass spectrometry. The multivariate nature of BMA combined with its ability to account for model uncertainty makes it an attractive candidate to extract predictive genes from any highdimensional biological data.
All analyses in this study were conducted using R statistical software http://www.rproject.org/. The Bioconductor packages for the iterative BMA algorithms for classification and survival analysis described in this paper are available for download from Bioconductor's website http://www.bioconductor.org/ as the iterativeBMA and iterativeBMAsurv packages respectively. Please visit our supplemental website for access to the breast cancer and DLBCL datasets, along with other helpful links and information.
Supplemental website
Software
Our software implementation is publicly available as a bioconductor package called "iterativeBMAsurv" http://www.bioconductor.org/packages/2.3/bioc/html/iterativeBMAsurv.html
Abbreviations
 AFT:

accelerated failure time
 BMA:

Bayesian Model Averaging
 BSS/WSS:

betweengroups to withingroups sum of squares ratio
 DLBCL:

diffuse large Bcell lymphoma
 iterativeBMAsurv:

iterative Bayesian Model Averaging for survival analysis
 LASSO:

least absolute shrinkage and selection operator
 PAM:

Prediction Analysis for Microarrays
 PLS:

partial least squares
 PPS:

partial predictive score
 RFE:

recursive feature elimination
 SAM:

Significance Analysis for Microarrays
 TNoM:

threshold number of misclassification score
Declarations
Acknowledgements
We would like to thank Drs. Isabelle Bichindaritz, Donald Chinn, Steve Hanks, Ian Painter, Deanna Petrochilos, and Chris Volinsky. Bumgarner is funded by NIHNHLBI P50 HL073996, NIHNIAID U54 AI057141, NIHNCRR R24 RR02186301A1, NIHNIDCR R01 DE01221206, NIHNCRR 1 UL1 RR 02501401, and a generous basic research grant from Merck. Raftery is supported by NIHNICHD 1R01HDO5451101A1, NSF llS0534094, NSF ATM0724721, and Office of Naval Research grant N000140110745. Yeung is supported by NIHNCI K25CA106988 and NIHNIGMS R01GM08416301A1.
Authors’ Affiliations
References
 Li J, Duan Y, Ruan X: A Novel Hybrid Approach to Selecting Marker Genes for Cancer Classification Using Gene Expression Data. The 1st International Conference on Bioinformatics and Biomedical Engineering, 2007, ICBBE. 2007, 264267.Google Scholar
 Liu H, Motoda H: Feature Selection for Knowledge Discovery and Data Mining. 1998, Boston: Kluwer Academic PublishersView ArticleGoogle Scholar
 Liu H, Motoda H: Computational Methods of Feature Selection. Chapman & Hall/CRC data mining and knowledge discovery series. 2008, Boca Raton: Chapman & Hall/CRC PressGoogle Scholar
 Nguyen D, Rocke D: Tumor classification by Partial Least Squares Using Microarray Gene Expression Data. Bioinformatics. 2002, 18: 3950. 10.1093/bioinformatics/18.1.39.View ArticlePubMedGoogle Scholar
 Golub T, Slonim D, Tamayo P, Huard C, Gaasenbeek M, Mesirov J, Coller H, Loh M, Downing J, Caliqiuri M, Bloomfield C, Lander E: Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring. Science. 1999, 286: 531537. 10.1126/science.286.5439.531.View ArticlePubMedGoogle Scholar
 Cox D: Regression Models and Life Tables. Journal of the Royal Statistical Society, Series B. 1972, 34: 187220.Google Scholar
 BenDor A, Bruhn L, Friedman N, Nachman I, Schummer M, Yakhini Z: Tissue Classification with Gene Expression Profiles. Journal of Computational Biology. 2000, 7: 559583. 10.1089/106652700750050943.View ArticlePubMedGoogle Scholar
 Dudoit S, Fridlyan J, Speed T: Comparison of Discrimination Methods for the Classification of Tumors Using Gene Expression Data. Journal of the American Statistical Association. 2002, 97: 7787. 10.1198/016214502753479248.View ArticleGoogle Scholar
 Chow M, Moler E, Mian I: Identifying Marker Genes in Transcription Profiling Data Using a Mixture of Feature Relevance Experts. Physiol Genomics. 2001, 5: 99111.PubMedGoogle Scholar
 Guyon I, Weston J, Barnhill S: Gene Selection for Cancer Classification Using Support Vector Machines. Machine Learning. 2002, 46: 389422. 10.1023/A:1012487302797.View ArticleGoogle Scholar
 Li L, Weinberg C, Darden T, Pedersen L: Gene Selection for Sample Classification Based on Gene Expression Data: Study of Sensitivity to Choice of Parameters of the GA/KNN Method. Bioinformatics. 2001, 17: 11311142. 10.1093/bioinformatics/17.12.1131.View ArticlePubMedGoogle Scholar
 Silva P, Hashimoto R, Kim S, Barrera J, Brandao L, Suh E, Dougherty E: Feature Selection Algorithms to Find Strong Genes. Pattern Recognition Letters. 2005, 26: 14441453. 10.1016/j.patrec.2004.11.017.View ArticleGoogle Scholar
 Yu J, Almal A, Dhanasekaran S, Ghosh D, Worzel W, Chinnaiyan A: Feature Selection and Molecular Classification of Cancer Using Genetic Programming. Neoplasia. 2007, 9: 292303. 10.1593/neo.07121.PubMed CentralView ArticlePubMedGoogle Scholar
 Pudil P, Novovicova J, Kittler J: Floating Search Methods in Feature Selection. Physical Review Letters. 1994, 15: 11191125.Google Scholar
 Geman D, D'Avignon C, Naiman D, Winslow R: Classifying Gene Expression Profiles from Pairwise mRNA Comparisons. Statistical Applications in Genetics and Molecular Biology. 2004, 3: 121. 10.2202/15446115.1071.View ArticleGoogle Scholar
 Xu L, Tan A, Naiman D, Geman D, Winslow R: Robust Prostate Cancer Marker Genes Emerge from Direct Integration of InterStudy Microarray Data. Bioinformatics. 2005, 21: 39053911. 10.1093/bioinformatics/bti647.View ArticlePubMedGoogle Scholar
 Lai C, Reinders M, van't Veer L, Wessels L: A Comparison of Univariate and Multivariate Gene Selection Techniques for Classification of Cancer Datasets. BMC Bioinformatics. 2006, 7: 23510.1186/147121057235.PubMed CentralView ArticlePubMedGoogle Scholar
 Chen C, Wu T, Wu Y, Huang Y, Lee J: Characterization of the Univariate and Multivariate Techniques on the Analysis of Simulated and fMRI Datasets with Visual Task. Nuclear Science Symposium Conference Record, 2003 IEEE. 2003, 4: 24682472.Google Scholar
 Yeung K, Bumgarner R, Raftery AE: Bayesian Model Averaging: Development of an Improved MultiClass, Gene Selection and Classification Tool for Microarray Data. Bioinformatics. 2005, 21: 23942402. 10.1093/bioinformatics/bti319.View ArticlePubMedGoogle Scholar
 Huang T, Kecman V, Kopriva I: Kernel Based Algorithms for Mining Huge Data Sets: Supervised, SemiSupervised, and Unsupervised Learning. Studies in computational intelligence. 2006, Berlin: Springer Verlag, 17:View ArticleGoogle Scholar
 Witten I, Frank R: Data mining: Practical Machine Learning Tools and Techniques. 2005, San Francisco: Elsevier, Inc, SecondGoogle Scholar
 Kotsiantis S: Supervised Machine Learning: A Review of Classification Techniques. Informatica. 2007, 31: 249268.Google Scholar
 Xu L, Geman D, Winslow R: LargeScale Integration of Cancer Microarray Data Identifies a Robust Common Cancer Signature. BMC Bioinformatics. 2007, 8: 27510.1186/147121058275.PubMed CentralView ArticlePubMedGoogle Scholar
 Jiang H, Deng Y, Chen H, Tao L, Sha Q, Chen J, Tsai C, Zhang S: Joint Analysis of Two Microarray GeneExpression Data Sets to Select Lung Adenocarcinoma Marker Genes. BMC Bioinformatics. 2004, 5: 8110.1186/14712105581.PubMed CentralView ArticlePubMedGoogle Scholar
 Tan A, Naiman D, Xu L, Winslow R, Geman D: Simple Decision Rules for Classifying Human Cancers from Gene Expression Profiles. Bioinformatics. 2005, 21: 38963904. 10.1093/bioinformatics/bti631.PubMed CentralView ArticlePubMedGoogle Scholar
 Sotiriou C, Neo S, McShane L, Korn E, Long P, Jazaeri A, Martiat P, Fox S, Harris A, Liu E: Breast Cancer Classification and Prognosis Based on Gene Expression Profiles from a PopulationBased Study. PNAS. 2003, 100: 1039310398. 10.1073/pnas.1732912100.PubMed CentralView ArticlePubMedGoogle Scholar
 van 't Veer LJ, Dai H, Vijver van de MJ, He YD, Hart AA, Mao M, Peterse HL, Kooy van der K, Marton MJ, Witteveen AT, Schreiber GJ, Kerkhoven RM, Roberts C, Linsley PS, Bernards R, Friend SH: Gene expression profiling predicts clinical outcome of breast cancer. Nature. 2002, 415: 530536. 10.1038/415530a.View ArticlePubMedGoogle Scholar
 Raponi M, Zhang Y, Yu J, Chen G, Lee G, Taylor J, MacDonald J, Thomas D, Moskaluk C, Wang Y, Beer D: Gene Expression Signatures for Predicting Prognosis of Squamous Cell and Adenocarcinomas of the Lung. Cancer Research. 2006, 66: 74667472. 10.1158/00085472.CAN061191.View ArticlePubMedGoogle Scholar
 Hu H, Li J, Plank A, Wang H, Daggard G: Comparative Study of Classification Methods for Microarray Data Analysis. Proceedings of the Fifth Australasian Conference on Data Mining and Analystics: 2006; Sydney, Australia. 2006, Australian Computer Society, Inc, 3337.Google Scholar
 Langley P: Selection of relevant features in machine learning. Proceedings of the AAAI Fall symposium on relevance: 1994. 1994, New Orleans: AAAI Press, 140144.Google Scholar
 Volinsky C, Madigan D, Raftery AE, Kronmal R: Bayesian Model Averaging in Proprtional Hazard Models: Assessing the Risk of a Stroke. Applied Statistics. 1997, 46: 443448.Google Scholar
 Hosmer DW, Lemeshow S, May S: Applied Survival Analysis: Regression Modeling of Time to Event Data. 2008, New York: John WileyView ArticleGoogle Scholar
 Beer D, Kardia S, Huang C, Giordano T, Levin A, Misek D, Lin L, Chen G, Gharib T, Thomas D, Lizyness M, Kuick R, Hayasaka S, Taylor J, Iannettoni M, Orringer M, Hanash S: GeneExpression Profiles Predict Survival of Patients with Lung Adenocarcinoma. Nature Medicine. 2002, 8: 816824.PubMedGoogle Scholar
 Lu Y, Lemon W, Liu P, Yi Y, Morrison C, Yang P, Sun Z, Szoke J, Gerald W, Watson M, Govindan R, You M: A Gene Expression Signature Predicts Survival of Patients with Stage I NonSmall Cell Lung Cancer. PLOS Medicine. 2006, 3: 22292243. 10.1371/journal.pmed.0030467.View ArticleGoogle Scholar
 Bair E, Tibshirani R: SemiSupervised Methods to Predict Patient Survival from Gene Expression Data. PLOS Biology. 2004, 2: 511522. 10.1371/journal.pbio.0020108.View ArticleGoogle Scholar
 Nguyen DV, Rocke DM: Partial least squares proportional hazard regression for application to DNA microarray survival data. Bioinformatics. 2002, 18: 16251632. 10.1093/bioinformatics/18.12.1625.View ArticlePubMedGoogle Scholar
 Li H, Gui J: Partial Cox regression analysis for highdimensional microarray gene expression data. Bioinformatics. 2004, 20 (Suppl 1): i208215. 10.1093/bioinformatics/bth900.View ArticlePubMedGoogle Scholar
 Tibshirani R: Regression shrinkage and selection via the LASSO. Journal of the Royal Statistical Society, Series B. 1996, 58: 267288.Google Scholar
 Tibshirani R: The lasso method for variable selection in the Cox model. Stat Med. 1997, 16: 385395. 10.1002/(SICI)10970258(19970228)16:4<385::AIDSIM380>3.0.CO;23.View ArticlePubMedGoogle Scholar
 Zhang HH, Lu W: Adaptive Lasso for Cox's proportional hazards model. Biometrika. 2007, 94: 691793. 10.1093/biomet/asm037.View ArticleGoogle Scholar
 Kaderali L, Zander T, Faigle U, Wolf J, Schultze JL, Schrader R: CASPAR: a hierarchical bayesian approach to predict survival times in cancer from gene expression data. Bioinformatics. 2006, 22: 14951502. 10.1093/bioinformatics/btl103.View ArticlePubMedGoogle Scholar
 Bovelstad HM, Nygard S, Storvold HL, Aldrin M, Borgan O, Frigessi A, Lingjaerde OC: Predicting survival from microarray data – a comparative study. Bioinformatics. 2007, 23: 20802087. 10.1093/bioinformatics/btm305.View ArticlePubMedGoogle Scholar
 Kalbfleisch JD, Prentice RL: The statistical analysis of failure time data. 1980, New York: WileyGoogle Scholar
 Huang J, Ma S, Xie H: Regularized estimation in the accelerated failure time model with highdimensional covariates. Biometrics. 2006, 62: 813820. 10.1111/j.15410420.2006.00562.x.View ArticlePubMedGoogle Scholar
 Cai T, Huang J, Tian L: Regularized Estimation for the Accelerated Failure Time Model. Biometrics. 2008,Google Scholar
 Datta S, LeRademacher J: Predicting patient survival from microarray data by accelerated failure time modeling using partial least squares and LASSO. Biometrics. 2007, 63: 259271. 10.1111/j.15410420.2006.00660.x.View ArticlePubMedGoogle Scholar
 Raftery AE: Bayesian Model Selection in Social Research (with Discussion). Sociological Methodology 1995. Edited by: Marsden P. 1995, Cambridge, MA: Blackwell, 111196. 10.2307/271063.Google Scholar
 Hoeting J, Madigan D, Raftery AE, Volinsky C: Bayesian Model Averaging: A Tutorial. Statistical Science. 1999, 14: 382417. 10.1214/ss/1009212519.View ArticleGoogle Scholar
 Vijver van de MJ, He YD, van't Veer LJ, Dai H, Hart AA, Voskuil DW, Schreiber GJ, Peterse JL, Roberts C, Marton MJ, Parrish M, Atsma D, Witteveen A, Glas A, Delahaye L, Velde van der T, Bartelink H, Rodenhuis S, Rutgers ET, Friend SH, Bernards R: A geneexpression signature as a predictor of survival in breast cancer. N Engl J Med. 2002, 347: 19992009. 10.1056/NEJMoa021967.View ArticlePubMedGoogle Scholar
 Rosenwald A, Wright G, Wing C, Connors J, Campo E, Fisher R, Gascoyne R, MullerHermelink H, Smeland E, Giltnane J, Hurt E, Zhao H, Averett L, Yang L, Wilson W, Jaffe E, Simon R, Klausner R, Powell J, Duffey P, Longo D, Greiner T, Weisenburger D, Sanger W, Dave B, Lynch J, Vose J, Armitage J, Montserrat E, LopezGuillermo A: The Use of Molecular Profiling to Predict Survival After Chemotherapy for Diffuse LargeBCell Lymphoma. The New England Journal of Medicine. 2002, 346: 19371947. 10.1056/NEJMoa012914.View ArticlePubMedGoogle Scholar
 A clinical evaluation of the International Lymphoma Study Group classification of nonHodgkin's lymphoma. The NonHodgkin's Lymphoma Classification Project. Blood. 1997, 89: 39093918.
 Shipp M, Ross K, Tamayo P, Weng A, Kutok J, Aguiar R, Gaasenbeek M, Angelo M, Reich M, Pinkus G, Ray T, Koval M, Last K, Norton A, Lister T, Mesirov J, Neuberg D, Lander E, Aster J, Golub T: Diffuse Large BCell Lymphoma Outcome Prediction by GeneExpression Profiling and Supervised Machine Learning. Nature Medicine. 2002, 8: 6874. 10.1038/nm010268.View ArticlePubMedGoogle Scholar
 Alizadeh A, Eisen M, Davis R, Ma C, Sabet H, Tran T, Powell J, Yang L, Marti G, Moore D, Hudson J, Chan W, Greiner T, Weisenburger D, Armitage J, Lossos I, Levy R, Botstein D, Brown P, Staudt L: The Lymphochip: A Specialized cDNA Microarray for the GenomicScale Analysis of Gene Expression in Normal and Malignant Lymphocytes. Cold Spring Harbor Symposia on Quantitative Biology. 1999, 64: 7178. 10.1101/sqb.1999.64.71.View ArticlePubMedGoogle Scholar
 Derksen S, Keselman H: Backward, Forward and Stepwise Automated Subset Selection Algorithms: Frequency of Obtaining Authentic and Noise Variables. British Journal of Mathematical and Statistical Psychology. 1992, 45: 265282.View ArticleGoogle Scholar
 Furnival G, Wilson R: Regression by Leaps and Bounds. Technometrics. 1974, 16: 499511. 10.2307/1267601.View ArticleGoogle Scholar
 Madigan D, Raftery AE: Model Selection and Accounting for Model Uncertainty in Graphical Models Using Occamís Window. Journal of the American Statistical Association. 1994, 89: 13351346. 10.2307/2291017.View ArticleGoogle Scholar
 Kuo L, Smith A: Bayesian Computations in Survival Models Via the Gibbs Sampler. Survival Analysis: State of the Art. Edited by: Klein J, Goel P. 1992, Boston: Dordrecht, 1124.View ArticleGoogle Scholar
 Draper D: Assessment and Propagation of Model Uncertainty. Journal of the Royal Statistical Society, Series B. 1995, 57: 4597.Google Scholar
 Taplin R: Robust Likelihood Calculation for Time Series. Journal of the Royal Statistical Society, Series B. 1993, 55: 829836.Google Scholar
 Taplin R, Raftery AE: Analysis of Agricultural Field Trials in the Presence of Outliers and Fertility Jumps. Biometrics. 1994, 50: 764781. 10.2307/2532790.View ArticleGoogle Scholar
 Volinsky C, Raftery AE: Bayesian Information Criterion for Censored Survival Models. Biometrics. 2000, 56: 256262. 10.1111/j.0006341X.2000.00256.x.View ArticlePubMedGoogle Scholar
 Raftery AE: Approximate Bayes Factors and Accounting for Model Uncertainty in Generalised Linear Models. Biometrika. 1996, 83: 251266. 10.1093/biomet/83.2.251.View ArticleGoogle Scholar
 Kaplan E, Meier P: Nonparametric Estimation from Incomplete Observations. Journal of the American Statistical Association. 1958, 53: 457481. 10.2307/2281868.View ArticleGoogle Scholar
 Supplementary web site for Predicting survival from microarray data – a comparative study. [http://www.med.uio.no/imb/stat/bmms/software/microsurv/]
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.