BACKGROUND
prediction methods  are widely used in biomedical research.
for example, reliable prediction methods are essential for accurate disease classification, diagnosis and prognosis.
since prediction methods based on multiple features  can greatly outperform predictors based on a single feature  <cit> , it is important to develop methods that can optimally combine features to obtain high accuracy.
introductory text books describe well known prediction methods such as linear discriminant analysis , k-nearest neighbor  predictors, support vector machines   <cit> , and tree predictors  <cit> .
many publications have evaluated popular prediction methods in the context of gene expression data  <cit> .
ensemble predictors are particularly attractive since they are known to lead to highly accurate predictions.
an ensemble predictor generates and integrates multiple versions of a single predictor , and arrives at a final prediction by aggregating the predictions of multiple base learners, e.g.
via plurality voting across the ensemble.
one particular approach for constructing an ensemble predictor is bootstrap aggregation   <cit> .
here multiple versions of the original data are generated through bootstrapping, where observations from the training set are randomly sampled with replacement.
an individual predictor  is fitted on each bootstrapped data set.
thus,  <dig> bootstrapped data sets  will lead to an ensemble of  <dig> tree predictors.
in case of a class outcome , the individual predictors “vote” for each class and the final prediction is obtained by majority voting.
breiman  showed that bagging weak predictors  often yields substantial gains in predictive accuracy  <cit> .
but it seems that ensemble predictors are only very rarely used for predicting clinical outcomes.
this fact points to a major weakness of ensemble predictors: they typically lead to ”black box” predictions that are hard to interpret in terms of the underlying features.
clinicians and epidemiologists prefer forward selected regression models since the resulting predictors are highly interpretable: a linear combination of relatively few features can be used to predict the outcome or the probability of an outcome.
but the sparsity afforded by forward feature selection comes at an unacceptably high cost: forward variable selection  often greatly overfit the data which results in unstable and inaccurate predictors  <cit> .
ideally, one would want to combine the advantages of ensemble predictors with those of forward selected regression models.
as discussed below, multiple articles describe ensemble predictors based on linear models including the seminal work by breiman  <cit>  who evaluated a bagged forward selected linear regression model.
however, the idea of bagging forward selected linear models  appears to have been set aside as new ensemble predictors, such as the random forest, became popular.
a random forest  predictor not only bags tree predictors but also introduces an element of randomness by considering only a randomly selected subset of features at each node split  <cit> .
the number of randomly selected features, mtry, is the only parameter of the random forest predictor.
the random forest predictor has deservedly received a lot of attention for the following reasons: first, the bootstrap aggregation step allows one to use out-of-bag  samples to estimate the predictive accuracy.
the resulting oob estimate of the accuracy often obviates the need for cross-validation and other resampling techniques.
second, the rf predictor provides several measures of feature  importance.
several articles explore the use of these importance measures to select genes  <cit> .
third, it can be used to define a dissimilarity measure that can be used in clustering applications  <cit> .
fourth, and most importantly, the rf predictor has superior predictive accuracy.
it performs as well as alternatives in cancer gene expression data sets  <cit>  but it really stands out when applied to the uci machine learning benchmark data sets where it is as good as  many existing methods  <cit> .
while we confirm the truly outstanding predictive performance of the rf, the proposed rglm method turns out to be even more accurate than the rf .
breiman and others have pointed out that the black box predictions of the rf predictor can be difficult to interpret.
for this reason, we wanted to give bagged forward selected generalized linear regression models another careful look.
after exploring different approaches for injecting elements of randomness into the individual glm predictors, we arrived at a new ensemble predictor, referred to as random glm predictor, with an astonishing predictive performance.
an attractive aspect of the proposed rglm predictor is that it combines the advantages of the rf with that of a forward selected glm.
as the name generalized linear model indicates, it can be used for a general outcome such as a binary outcome, a multi-class outcome, a count outcome, and a quantitative outcome.
we show that several incremental  changes to the original bagged glm predictor by breiman add to up to a qualitatively new predictor  that performs at least as well as the rf predictor on the uci benchmark data sets.
while the uci data are the benchmark data for evaluating predictors, only a dozen such data sets are available for binary outcome prediction.
to provide a more comprehensive empirical comparison of the different prediction methods, we also consider over  <dig> comparisons involving gene expression data.
in these genomic studies, the rglm method turns out to be slightly more accurate than the considered alternatives.
while the improvements in accuracy afforded by the rglm are relatively small they are statistically significant.
this article is organized as follows.
first, we present a motivating example that illustrates the high prediction accuracy of the rglm.
second, we compare the rglm with other state of the art predictors when it comes to binary outcome prediction.
toward this end, we use the uci machine learning benchmark data, over  <dig> empirical gene expression comparisons, and extensive simulations.
third, we compare the rglm with other predictors for quantitative  outcome prediction.
fourth, we describe several variable importance measures and show how they can be used to define a thinned version of the rglm that only uses few important features.
even for data sets comprised of thousands of gene features, the thinned rglm often involves fewer than  <dig> features and is thus more interpretable than most ensemble predictors.
methods
construction of the rglm predictor
rglm is an ensemble predictor based on bootstrap aggregation  of generalized linear models whose features  are selected using forward regression according to aic criterion.
glms comprise a large class of regression models, e.g.
linear regression for a normally distributed outcome, logistic regression for binary outcome, multi-nomial regression for multi-class outcome and poisson regression for count outcome  <cit> .
thus, rglm can be used to predict binary-, continuous-, count-, and other outcomes for which generalized linear models can be defined.
the “randomness” in rglm stems results from two sources.
first, a non-parametric bootstrap procedure is used which randomly selects samples with replacement from the original data set.
second, a random subset of features  is selected for each bootstrap sample.
this amounts to a random sub-space method  <cit>  applied to each bootstrap sample separately.
the steps of the rglm construction are presented in figure  <dig>  first, starting from the original data set another equal-sized data set is generated using the non-parametric bootstrap method, i.e.
samples are selected with replacement from the original data set.
the parameter nbags  determines how many of such bootstrap data sets  are being generated.
second, a random set of features  is randomly chosen for each bag.
thus, the glm predictor for bag  <dig> will typically involve a different set of features than that for bag  <dig>  third, the nfeaturesinbag of randomly selected features per bag are rank-ordered according to their individual association with the outcome variable y in each bag.
for a quantitative outcome y, one can simply use the absolute value of the correlation coefficient between the outcome and each feature to rank the features.
more generally, one can fit a univariate glm model to each feature to arrive at an association measure .
only the top ranking features  will become candidate covariates for forward selection in a multivariate regression model.
the top number of candidate features is determined by the input parameter ncandidatecovariates .
fourth, forward variable  selection is applied to the ncandidatecovariates of each bag to arrive at a multivariable generalized linear model per bag.
the forward selection procedure used by rglm is based on the stepaic r function in the mass r library where method is set to “forward”.
fifth, the predictions of each forward selected multivariate model  are aggregated across bags to arrive at a final ensemble prediction.
the aggregation method depends on the type of outcome.
for a continuous outcome, predicted values are simply averaged across bags.
for a binary outcome, the adjusted majority vote  strategy  <cit>  is used which averages predicted probabilities across bags.
given the estimated class probabilities one can get a binary prediction by choosing an appropriate threshold .
importantly, rglm also has a parameter maxinteractionorder  for creating interactions up to a given order among features in the model construction.
for example, rglm.inter <dig> results from setting maxinteractionorder= <dig>  i.e.
considering pairwise  interaction terms.
as example, consider the case when only pairwise interaction terms are used.
for each bag a random set of features is selected  from the original covariates, i.e.
covariates without interaction terms.
next, all pairwise interactions among the nfeaturesinbag randomly selected features are generated.
next, the usual rglm candidate feature selection steps will be applied to the combined set of pairwise interaction terms and the nfeaturesinbag randomly selected features per bag resulting in ncandidatecovariates top ranking features per bag, which are subsequently subjected to forward feature selection.
these methods are implemented in our r software package randomglm which allows the user to input a training set and optionally a test set.
it automatically outputs out-of-bag estimates of the accuracy and variable importance measures.
parameter choices for the rglm predictor
as discussed below, we find that it is usually sufficient to consider only nbags= <dig> bootstrap data sets.
the default value for nfeaturesinbag depends on the total number of features.
it is easier to explain it in terms of the proportion of features randomly selected per bag, nfeaturesinbag/n, where n is the total number of features of the training set.
apart from n it is also valuable to consider the effective number of features which equals the number of features n plus the number of interaction terms, e.g.
n∗=n+n/ <dig> in case of pairwise interactions.
using this notation, the default value of nfeaturesinbag can be arrived at by solving equations presented in table  <dig>  these equations were found by empirically evaluating various choices of nfeaturesinbag values .
in particular, we found that in case of n∗<= <dig>  then using all features  is often a good choice, whereas if n∗> <dig> then setting nfeaturesinbag/n= <dig>  works well.
the default value nfeaturesinbag/n= <dig> − <dig> n∗ in the intermediate case  results from fitting an interpolation line through the two points  and .
we find that rglm is quite robust with respect to the parameter nfeaturesinbag.
to limit the number of covariates considered in forward selection , the default value of ncandidatecovariates is set to  <dig>  overall, the default values perform well in our simulations, empirical gene expression and machine learning benchmark studies.
but we recommend to use the oob estimate of predictive accuracy to inform the choice of the parameter values.
this table shows the default values of nfeaturesinbag in terms of nfeaturesinbag/n for rglm, rglm.inter <dig> and rglm.inter <dig>  n is the total number of features of the training data.
n∗ is the effective number of features which equals the number of features n plus the number of interaction terms.
formulas are shown in terms of both n and the corresponding n∗.
<dig>  and  <dig>  are obtained by interpolating a straight line between  and .
relationship with related prediction methods
as discussed below, rglm can be interpreted as a variant of a bagged predictor  <cit> .
in particular, it is similar to the bagged forward linear regression model  <cit>  but differs in the following aspects:
<dig>  rglm allows for interaction terms between features which greatly improve the performance on some data sets .
we refer to rglm involving two-way or three way interactions as rglm.inter <dig> and rglm.inter <dig>  respectively.
<dig>  rglm has a parameter nfeaturesinbag that allows one to restrict the number of features used in each bootstrap sample.
this parameter is conceptually related to the mtry parameter of the random forest predictor.
in essence, this parameter allows one to use a random subspace method  in each bootstrap sample.
<dig>  rglm has a parameter ncandidatecovariates that allows one to restrict the number of features in forward regression, which not only has computational advantages but also introduces additional instability into the individual predictors, which is a desirable characteristic of an ensemble predictor.
<dig>  rglm optimizes the aic criterion during forward selection.
<dig>  rglm has a “thinning threshold” parameter which allows one to reduce the number of features involved in prediction while maintaining good prediction accuracy.
since a thinned rglm involves far fewer features, it facilitates the understanding how the ensemble arrives at its predictions.
rglm is not only related to bagging but also to the random subspace method  proposed by  <cit> .
in the rsm, the training set is also repeatedly modified as in bagging but this modification is performed in the feature space .
in the rsm, a subset of features is randomly selected which amounts to restricting attention to a subspace of the original feature space.
as one of its construction steps, rglm uses a rsm on each bootstrap sample.
future research could explore whether random partitions as opposed to random subspaces would be useful for constructing an rglm.
random partitions of the feature space are similar to random subspaces but they divide the feature space into mutually exclusive subspaces  <cit> .
random partition based predictors have been shown to perform well in high-dimensional data .
both rsm and random partitions have more general applicability than rglm since these methods can be used for any base learner.
there is a vast literature on ensemble induction methods but a property worth highlighting is that rglm uses forward variable selection of glms.
recall that rglm goes through the following steps: 1) bootstrap sampling, 2) rsm , 3) forward variable selection of a glm, 4) aggregation of votes.
empirical studies involving different base learners  have shown that combining bootstrap sampling with rsm  leads to ensemble predictors with comparable performance to that of the random forest predictor  <cit> .
another prediction method, random multinomial logit model , also shares a similar idea with rglm.
it was recently proposed for multi-class outcome prediction  <cit> .
rmnl bags multinomial logit models with random feature selection in each bag.
it can be seen as a special case of rglm, except that no forward model selection is carried out.
software implementation
the rglm method is implemented in the freely available r package randomglm.
the r function randomglm allows the user to output training set predictions, out-of-bag predictions, test set predictions, coefficient values, and variable importance measures.
the predict function can be used arrive at test set predictions.
tutorials can be found at the following webpage: http://labs.genetics.ucla.edu/horvath/rglm.
short description of alternative prediction methods
forward selected generalized linear model predictor  we denote by forwardglm the  generalized linear model predictor whose covariates were selected using forward feature selection .
thus, forwardglm does not involve bagging, random feature selection, and is not an ensemble predictor.
random forest  rf is an ensemble predictor that consists of a collection of decision trees which vote for the class of observations  <cit> .
the rf is known for its outstanding predictive accuracy.
we used the randomforest r package in our studies.
we considered two choices for the rf parameter mtry: i) the default rf predictor where mtry equals the square root of the number of features and ii) rfbigmtry where mtry equals the total number of features.
we always generated at least  <dig> trees per forest but used  <dig> trees when calculating variable importance measures.
recursive partitioning and regression trees  classification and regression trees were generated using the default settings rpart r package.
tree methods are described in  <cit> .
linear discriminant analysis  lda aims to find a linear combination of features  to predict a binary outcome .
we used the lda r function in the mass r package with parameter choice method=moment.
diagonal linear discriminant analysis  dlda is similar to lda but it ignores the correlation patterns between features.
while this is often an unrealistic assumption, dlda  has been found to work well in in gene expression applications  <cit> .
here we used the default parameters from the supclust r package  <cit> .
k nearest neighbor  we used the knn r function in the class r package  <cit> , which chose the parameter k of nearest neighbors using 3-fold cross validation .
support vector machines  we used the default parameters from the e <dig> r package to fit svms  <cit> .
additional details can be found in  <cit> .
shrunken centroids  the sc predictor is known to work well in the context of gene expression data  <cit> .
here we used the implementation in the pamr r package  <cit>  which chose the optimal level of shrinkage using cross validation.
penalized regression models various convex penalties can be applied to generalized linear models.
we considered ridge regression  <cit>  corresponding to an ℓ <dig> penalty, the lasso corresponding to an ℓ <dig> penalty  <cit> , and elastic net corresponding to a linear combination of ℓ <dig> and ℓ <dig> penalties  <cit> .
we used the glmnet r function from the glmnet r package  <cit>  with alpha parameter values of  <dig>   <dig>  and  <dig>  respectively.
glmnet also involves another parameter  which was chosen as the median of the lambda sequence output resulting from glmnet.
for uci benchmark data sets, pairwise interaction between features were considered.
<dig> disease-related gene expression data sets
we use  <dig> disease related gene expression data sets involving cancer and other human diseases .
the first  <dig> data sets involving various cancers were previously used by  <cit> .
these data can be downloaded from the author’s webpage at http://ligarto.org/rdiaz/papers/rfvs/randomforestvarsel.html.
the braintumor <dig> and dlbcl data sets were downloaded from http://www.gems-system.org/.
the remaining  <dig> data sets  were downloaded from either the gene expression omnibus  database or the arrayexpress data base in raw form and subsequently preprocessed using mas <dig> normalization and quantile normalization.
only the top  <dig> probes  with highest mean expression were considered for outcome prediction.
we briefly point out that diaz et al  report prediction error rates estimated using a bootstrap method.
in contrast, we report 3-fold cross validation estimates , which may explain minor numerical differences between our study and that of diaz et al .
sample size, number of features, original reference, data set ids and outcomes for the  <dig> disease related gene expression data sets.
empirical gene expression data sets
for all data sets below, we considered  <dig> randomly selected gene traits, i.e.
<dig> randomly selected probes.
they were directly used as continuous outcomes or dichotomized according to the median value  to generate binary outcomes.
for all data sets except “brain cancer”,  <dig> of the observations  were randomly chosen as the training set, while the remaining samples were chosen as test set.
we focused on the  <dig> genes  with the highest mean expression levels in each data set.
brain cancer data sets these two related data sets contain  <dig> and  <dig> microarray samples of glioblastoma  patients, respectively.
gene expression profiles were measured using affymetrix u <dig> microarrays.
a detailed description can be found in  <cit> .
the first data set  was used as a training set while and the second data set  was used as a test set.
safhs blood lymphocyte data set this data set  <cit>  was derived from blood lymphocytes of randomly ascertained participants enrolled in the san antonio family heart study.
gene expression profiles were measured with the illumina sentrix human whole genome  series i beadchips.
after removing potential outliers ,  <dig> samples remained in the data set.
wb whole blood gene expression data set this is the whole blood gene expression data from healthy controls.
peripheral blood samples from healthy individuals were analyzed using illumina human ht- <dig> microarrays.
after pre-processing,  <dig> samples remained in the data set.
mouse tissue gene expression data sets the  <dig> tissue specific gene expression data sets were generated by the lab of jake lusis at ucla.
these data sets measure gene expression levels  from adipose , brain , liver  and muscle  tissue of mice from the b ×h f <dig> mouse intercross described in  <cit> .
in addition to gene traits, we also predicted  <dig> quantitative mouse clinical traits including mouse weight, length, abdominal fat, other fat, total fat, adiposity index , plasma triglycerides, total plasma cholesterol, high-density lipoprotein fraction of cholesterol, plasma unesterified cholesterol, plasma free fatty acids, plasma glucose, plasma low-density lipoprotein and very low-density lipoprotein cholesterol, plasma mcp- <dig> protein levels, plasma insulin, plasma glucose-insulin ratio, plasma leptin, plasma adiponectin, aortic lesion size , aneurysms , and aortic calcification in the lesion area.
machine learning benchmark data sets
the  <dig> machine learning benchmark data sets used in this article are listed in table  <dig>  note that only eight of the  <dig> data sets have a binary outcomes.
the multi-class outcomes of the  <dig> remaining data sets were turned into binary outcomes by considering the most prevalent class versus all other classes combined.
missing data were imputed using nearest neighbor averaging.
for each data set and prediction method, we report the average 3-fold cv estimate of prediction accuracy over  <dig> random partitions of the data into  <dig> folds.
sample size and number of features for the  <dig> uci machine learning benchmark data sets.
simulated gene expression data sets
we simulated an outcome variable y and gene expression data that contained  <dig> modules .
only  <dig> of the modules were comprised of genes that correlated with the outcome y.
45% of the genes were background genes, i.e.
these genes were outside of any module.
the simulation scheme is detailed in additional file  <dig> and implemented in the r function simulatedatexpr5modules from the wgcna r package  <cit> .
this r function was used to simulate pairs of training and test data sets.
the simulation study was used to evaluate prediction methods for continuous outcomes and for binary outcomes.
for binary outcome prediction, the continuous outcome y was thresholded according to its median value.
we considered  <dig> different simulation scenarios involving varying sizes of the training data  and varying numbers of genes  that served as features.
test sets contained the same number of genes as in the corresponding training set and  <dig> samples.
for each simulation scenario, we simulate  <dig> replicates resulting from different choices of the random seed.
RESULTS
motivating example: disease-related gene expression data sets
we compare the prediction accuracy of rglm with that of other widely used methods on  <dig> gene expression data sets involving human disease related outcomes.
many of the  <dig> data sets  are well known cancer data sets, which have been used in other comparative studies  <cit> .
a brief description of the data sets can be found in methods.
to arrive at an unbiased estimate of prediction accuracy, we used 3-fold cross validation .
note that the accuracy equals  <dig> minus the median misclassification error rate.
table  <dig> reports the prediction accuracy of different methods including rglm, random forest , random forest , tree predictor , linear discriminant analysis , diagonal linear discriminant analysis , k nearest neighbor , support vector machine  and shrunken centroid .
a short description of these prediction methods is provided in methods.
for each data set, the prediction accuracy was estimated using 3−fold cross validation across  <dig> random partitions of the data into  <dig> folds.
mean accuracies across the  <dig> data sets and the resulting ranks are summarized at the bottom.
two sided paired wilcoxon test p-values can be used to determine whether the accuracy of rglm is significantly different from that of other predictors.
note that the rglm yields the highest mean accuracy.
as seen from table  <dig>  rglm achieves the highest mean accuracy in these disease data sets, followed by rfbigmtry and sc.
note that the standard random forest predictor  performs worse than rglm.
the accuracy difference between rglm and alternative methods is statistically significant  for all predictors except for rfbigmtry, dlda and sc.
since rfbigmtry is an ensemble predictor that relies on thousands of features it would be difficult to interpret its predictions in terms of the underlying genes.
our evaluations focused on the accuracy .
however, a host of other accuracy measures could be considered.
additional file  <dig> presents the results for sensitivity and specificity.
the top  <dig> methods with highest sensitivity are: rf , svm  and rglm .
the top  <dig> methods with highest specificity are: sc , rglm  and knn .
a strength of this empirical comparison is that it involves clinically or biologically interesting data sets but a severe limitation is that it only involves  <dig> comparisons.
therefore, we now turn to more comprehensive empirical comparisons.
binary outcome prediction
empirical study involving dichotomized gene traits
many previous empirical comparisons of gene expression data considered fewer than  <dig> data sets.
to arrive at  <dig> comparisons, we use the following approach: we started out with  <dig> human and mouse gene expression data sets.
for each data set, we randomly chose  <dig> genes as gene traits  resulting in 7× <dig> possible outcomes.
we removed the gene corresponding to the gene trait from the feature set.
next, each gene trait was dichotomized by its median value to arrive at a binary outcome y. the goal of each prediction analysis was to predict the dichotomized gene trait y based on the other genes.
at first sight, this artificial outcome is clinically uninteresting but it is worth emphasizing that clinicians often deal with dichotomized measures of gene products, e.g.
high serum creatinine levels may indicate kidney damage, high psa levels may indicate prostate cancer, and high hdl levels may indicate hypercholesterolemia.
to arrive at unbiased estimates of prediction accuracy, we split each data set into a training and test set.
figure  <dig>  shows boxplots of the accuracies across the  <dig> comparisons.
similar performance patterns are observed for the individual data sets ).
the figure also reports pairwise comparisons of the rglm method versus alternative methods.
specifically, it reports the two-sided wilcoxon signed rank test p-values for testing whether the accuracy of the rglm predictor is higher than that of the considered alternative method.
strikingly, rglm is more accurate than the other methods overall.
while the increase in accuracy is often minor, it is statistically significant as can be seen by comparing rglm to rf , rfbigmtry , lda , svm  and sc .
other predictors perform even worse, and the corresponding p-values are not shown.
the fact that rfbigmtry is more accurate in this situation than the default version of rf probably indicates that relatively few genes are informative for predicting a dichotomized gene trait.
also note that rglm is much more accurate than the unbagged forward selected glm which reflects that forward selection greatly overfits the training data.
in conclusion, these comprehensive gene expression studies show that rglm has outstanding prediction accuracy.
machine learning benchmark data analysis
here we evaluate the performance of rglm on the uci machine learning benchmark data sets which are often used for evaluating prediction methods  <cit> .
we consider  <dig> benchmark data sets from the mlbench r package:  <dig> uci data sets and  <dig> synthetic data sets .
we choose these data sets for two reasons.
first, these  <dig> data sets were also used in the original evaluation of the random forest predictor  <cit> .
second, these data include all of the available data sets with binary outcomes in the mlbench r package.
a detailed description of these data sets can be found in methods.
in his original publication on the random forest, breiman found that the rf outperformed bagged predictors on the uci benchmark data which may explain why bagged glms have not received much attention.
we hypothesize that the relatively poor performance of a bagged logistic regression model on these data sets could be ameliorated by considering interaction terms between the features.
table  <dig> confirms our hypothesis.
rglm.inter <dig>  has superior or tied accuracy compared to rglm in  <dig> out of  <dig> benchmark data sets.
in particular, pairwise interactions greatly improve the prediction accuracy in the ringnorm data set.
higher order interactions  do not perform better than rglm.inter <dig> but dramatically increase computational burden .
for each data set, the prediction accuracy was estimated using 3−fold cross validation across  <dig> random partitions of the data into  <dig> folds.
rglm.inter <dig> incorporates pairwise interaction between features into the rglm predictor.
mean accuracies and the resulting ranks are summarized at the bottom.
the wilcoxon signed rank test was used to test whether accuracy differences between rglm.inter <dig> and other predictors are significant.
rglm.inter <dig>  rf, and svm tie for first place .
overall, we find that rglm.inter <dig> ties with svm  and rf  for the first place in the benchmark data.
as can be seen from additional file  <dig>  rglm.inter <dig> achieves the highest sensitivity and specificity, which also support its good performance in the benchmark data sets.
a potential limitation of these comparisons is that we considered pairwise interaction terms for the rglm predictor but not for the other predictors.
to address this issue, we also considered pairwise interactions among features for other predictors.
additional file  <dig> shows that no method surpasses rglm.inter <dig> when pairwise interaction terms are considered.
in particular, interaction terms between features do not improve the performance of the random forest predictor.
a noteworthy disadvantage of rglm.inter in case of many features is the computational burden that may result from adding interaction terms.
in applications where interaction terms are needed for rglm, faster alternatives  remain an attractive choice.
simulation study involving binary outcomes
as described in methods, we simulated  <dig> gene expression data sets with binary outcomes.
the number of features  ranged from  <dig> to  <dig>  the sample sizes  of the training data ranged from  <dig> to  <dig>  to robustly estimate the test set accuracy we chose a large size for the corresponding test set data, n= <dig>  figure  <dig> shows the boxplots of the test set accuracies of different predictors.
the accuracy of the forwardglm is much lower than that of rglm, demonstrating the benefit of creating an ensemble predictor.
the wilcoxon test p-value shows that rglm is significantly better than all other methods except for the rf .
in this simulation study, rglm takes the first place when it comes to the median accuracy.
continuous outcome prediction
in the following, we show that rglm also performs exceptionally well when dealing with continuous quantitative outcomes.
we not only compare rglm to a standard forward selected linear model predictor  but also a random forest predictor .
we do not report the findings for the k-nearest neighbor predictor of a continuous outcome since it performed much worse than the above mentioned approaches in our gene expression applications .
we again split the data into training and test sets.
we use the correlation between test set predictions and truly observed test set outcomes as measure of predictive accuracy.
note that this correlation coefficient can take on negative values .
empirical study involving continuous gene traits
here we used the same  <dig> gene expression comparisons as described above  but did not dichotomize the gene traits.
incidentally, prediction methods for gene traits are often used for imputing missing gene expression values.
our results presented in figure  <dig> indicate that for the majority of genes high accuracies can be achieved.
but for some gene traits, the accuracy measure, which is defined as a correlation coefficient, takes on negative values indicating that there is no signal in the data.
note that the forward selected linear predictor ties with the random forest irrespective of the choice of the mtry parameter and both methods perform significantly worse than the rglm predictor.
mouse tissue expression data involving continuous clinical outcomes
here we used the mouse liver and adipose tissue gene expression data sets to predict  <dig> clinical outcomes .
again, rglm achieved significantly higher median prediction accuracy compared to the other predictors .
simulation study involving continuous outcomes
<dig> gene expression data sets are simulated in the same way as described previously  but here the outcome y was not dichotomized.
as shown in figure  <dig>  the forwardglm accuracy trails both rglm and rf, reflecting again the fact that forward regression overfits the data.
in this simulation study, we find that rglm yields significantly higher prediction accuracy than other predictors.
comparing rglm with penalized regression models
in our previous comparisons, we found that rglm greatly outperforms forward selected glm methods based on the aic criterion.
many powerful alternatives to forward variable selection have been developed in the literature, in particular penalized regression models.
here, we compare rglm to  <dig> major types of penalized regression models: ridge regression  <cit> , elastic net  <cit> , and the lasso  <cit> .
the predictive accuracies of these penalized regression models were compared to those of the rglm predictor using the same data sets described above for evaluating binary outcome and quantitative outcome prediction methods.
wilcoxon’s signed rank test was used to determine whether differences in predictive accuracy were significant.
figure  <dig>  shows that rglm outperforms penalized regression models when applied to binary outcomes.
for all comparisons, the paired median difference  is positive which indicates that rglm is at least as good if not better than any of these  <dig> penalized regression models.
in particular, rglm is significantly better than ridge regression  and the lasso  on the  <dig> dichotomized gene expression trait data.
also, rglm is significantly better than elastic net  and lasso  in simulations with binary outcomes.
figure  <dig>  shows that rglm outperforms penalized regression models for continuous outcome prediction as well.
positive accuracy differences again imply that rglm is at least as good as these penalized regression models.
in particular, it significantly outperforms ridge regression  in the  <dig> continuous gene expression traits data and outperforms elastic net  and lasso  in simulations with continuous outcomes.
as a caveat, we mention that cross validation methods were not used to inform the parameter choices of the penalized regression models since the rglm predictor was also not allowed to fine tune its parameters.
by only using default parameter choices we ensure a fair comparison.
in a secondary analysis, however, we allowed penalized regression models to use cross validation for informing the choice of the parameters.
while this slightly improved the performance of the penalized regression models , it did not affect our main conclusion.
rglm outperforms penalized regression models in these comparisons.
feature selection
here we briefly describe how rglm naturally gives rise to variable  importance measures.
we compare the variable importance measures of rglm with alternative approaches and show how variable importance measures can be used for defining a thinned rglm predictor with few features.
variable importance measure
there is a vast literature on using ensemble predictors and bagging for selecting features.
for example, meinshausen and bühlmann describe “stability selection” based on variable selection employed in regression models  <cit> .
the method involves repetitive sub-sampling, and variables that occur in a large fraction of the resulting selection set are chosen.
li et al.
use a random k-nearest neighbor predictor  to carry out feature selection  <cit> .
the entropy-based recursive feature elimination  method of furlanello et al.
ranks features in high dimensional microarray data  <cit> .
rglm, like many ensemble predictors, gives rise to several measures of feature  importance.
for example, the number of times a feature is selected in the forward glm across bags, timesselectedbyforwardregression, is a natural measure of variable importance .
another variable importance measure is the number of times a feature is selected as candidate covariate for forward regression, timesselectedascandidates.
note that both timesselectedbyforwardregression and timesselectedascandidates have to be ≤nbags.
finally, one can use the sum of absolute glm coefficient values, sumabscoefbyforwardregression, as a variable importance measure.
we prefer timesselectedbyforwardregression, since it is more intuitive and points to the features that directly contribute to outcome prediction.
to reveal relationships between different types of variable importance measures, we present a hierarchical cluster tree of rglm measures, rf measures and standard marginal analysis based on correlations in figure  <dig>  as expected, the marginal association measures  cluster together.
the same holds for the random forest based importance measures  and the  <dig> rglm based importance measures.
leo breiman already pointed out that random forests could be used for feature selection in genomic applications.
díaz-uriarte et al.
proposed a related gene selection method based on the rf which yields small sets of genes  <cit> .
this rf based gene selection method does not return sets of genes that are highly correlated because such genes would be redundant when it comes to predictive purposes.
since the rglm based importance measure timesselectedbyforwardregression is expected to lean towards selecting genes that are highly associated with the outcome, it comes as no surprise that only a few genes selected by the procedure of díaz-uriarte et al.
turn out to have a top ranking in terms of the rglm measure timesselectedbyforwardregression .
it is beyond our scope to provide a thorough evaluation of the different variable selection approaches and we refer the reader to the literature, e.g.
<cit> .
while our studies show that the rglm based variable importance measures have some relationships to other measures, they are sufficiently different from other measures to warrant a thorough evaluation in future comparison studies.
rglm predictor thinning based on a variable importance measure
both rglm and random forest have superior prediction accuracy but they differ with respect to how many features are being used.
recall that the random forest is composed of individual trees.
each tree is constructed by repeated node splits.
the number of features considered at each node split is determined by the rf parameter mtry.
the default value of mtry is the square root of the number of features.
in case of  <dig> gene features in our empirical studies, the default value is mtry= <dig>  for rfbigmtry, we choose all possible features, i.e.
mtry= <dig>  we find that a random forest predictor typically uses more than 40% of the features  in the empirical studies.
in contrast, rglm typically only involves a few hundred genes in these studies.
there are several reasons why rglm uses far fewer features in its construction.
first, and foremost, it uses forward selection  to select features in each bag.
second, the number of candidate covariates considered for forward regression is chosen to be low, i.e.
ncandidatecovariates= <dig>
in rglm, the number of times a feature is selected by forward regression models among all bags, timesselectedbyforwardregression, follows a highly skewed distribution.
only few features are repeatedly selected into the model while most features are selected only once .
it stands to reason that an even sparser, highly accurate predictor can be defined by refitting the glm on each bag without considering these rarely selected features.
we refer to this feature removal process as rglm predictor thinning.
thus, features whose value of timesselectedbyforwardregression lies below a pre-specified thinning threshold will be removed from the model fit a posteriori.
interestingly, the accuracy diminishes very slowly for initial, low threshold values.
but even low threshold values lead to a markedly sparser ensemble predictor ).
in other words, the average fraction of features  remaining in the thinned rglm declines drastically as the thinning threshold increases.
we have found that the following empirical function accurately describes the relationship between thinning threshold  and proportion of features left in the thinned rglm predictor:
proportionleft=f=1x=0exp{−e <dig> nbags <dig> )}0<x≤ <dig>
where x=thinningthresholdnbags and e denotes euler’s constant e≈ <dig> .
equation  <dig> was found by log transforming the data and using optimization approaches for estimating the parameters.
no mathematical derivation was used.
one can easily show that f  is a monotonically decreasing function which accurately describes the proportion of remaining features as can be seen from figure  <dig> .
since the proportion of remaining variables depends not only on the thinning threshold but also on the number of bags nbags, we also study how these results depend on the choice of nbags.
toward this end, we varied nbags from  <dig> to  <dig> for predicting the  <dig> dichotomized gene traits in the mouse adipose data set.
additional file  <dig> shows that the predicted values  based on equation  <dig> overlaps almost perfectly with the observed values  for all considered choices of nbags, which indicates that equation  <dig> accurately estimates the proportion of remaining features for range of different values of nbags.
our results demonstrate that the number of required features decreases rapidly even for low values of the thinning threshold without compromising the prediction accuracy of the thinned predictor.
figure  <dig>  shows that a thinning threshold of  <dig>  leads to a thinned predictor whose accuracy is negligibly lower  than that of the original rglm predictor but it involves less than 20% of the original number of variables.
recall that even the original number of variables is markedly lower than that of the rf predictor.
these results demonstrate that the thinned rglm combines the advantages of an ensemble predictor  with that of a forward selected glm model .
rglm thinning versus rf thinning
the idea behind rglm thinning is to remove features with low values of the variable importance measure.
of course, a similar idea can be applied to other predictors.
here we briefly evaluate the performance of a thinned random forest predictor which removed variables based on a low value of its importance measure .
to arrive at an unbiased comparison, both rglm and rf are thinned based on results obtained in the training data.
next, accuracies of the thinned predictors are evaluated in the test set data.
figure  <dig> compares thinned rglm versus thinned rf in our disease related data sets and also the empirical studies.
numbers that connect dashed lines are rglm thinning thresholds.
for a pre-specified threshold, the number of features used in the thinned random forest is matched to that used in the thinned rglm .
without thinning, rf uses a lot more features than rglm as mentioned previously.
as expected, the median number of genes left for prediction and the corresponding median prediction accuracy generally decrease as the thinning threshold becomes more stringent.
overall, a thinned rglm yields a significantly higher median accuracy than a thinned rf across different thinning thresholds .
in clinical practice, a thinned predictor with very few features and good accuracy can be very useful and interpretable.
for example, choosing a threshold of  <dig> in panel figure  <dig>  and a threshold of  <dig> in panel  would result in very sparse predictors.
in both cases, especially in panel , the thinned rglm has higher median accuracy than that of the thinned rf.
discussion
why was the rglm not discovered earlier?
after breiman proposed the idea of bagged linear regression models in  <dig>  <cit> , many authors have explored the utility of bagging logistic regression models  <cit> .
most previous studies report that bagging does not improve the accuracy of logistic regression.
bühlman and yu showed theoretically that bagging helps for “hard threshold” methods but not for “soft threshold” methods   <cit> .
these studies indicate that bagged logistic regression models are not beneficial since the individual predictors  are too stable.
overall, we agree with these results.
but our comprehensive evaluations show that by injecting elements of randomness and instability into a bagged logistic regression model one arrives at a state of the art prediction method that often outperforms existing methods.
figure  <dig> describes why the construction of the rglm runs counter to conventional wisdom.
as indicated by the upper right hand panel of figure  <dig>  the rglm is based on two seemingly bad modifications to a glm.
as indicated by the top left panel of figure  <dig>  forward selection of a glm is typically a bad idea since it overfits the data and thus degrades the prediction accuracy of a single glm predictor.
as indicated by the bottom right panel of figure  <dig>  bagging a full logistic regression  is also a bad idea since it leads to a complicated  predictor without clear evidence for increased accuracy .
but these two seemingly bad modifications add up to a superior prediction method .
breiman already noted that the instability afforded by variable selection is important for constructing a bagged linear model based predictor  <cit> .
in order to define an accurate glm based ensemble predictor, we also find that it is important to introduce additional elements of randomness and instability, which is also reflected in the name random glm.
our results show that the proposed changes  results in a more accurate predictor that involves surprisingly few features .
additional reasons why the merits of rglm have not been recognized earlier may be the following.
first, it may be a historical accident.
bagging was quickly over-shadowed by other seemingly more accurate ways of constructing ensemble predictors, such as boosting  <cit>  and the rf  <cit> , both of which have markedly better performance on the uci benchmark data.
we find that rglm.inter <dig> ties with svm and rf for the top spot in uci benchmark data set .
incidentally, rglm performs significantly better than svm and rf on the disease data sets  and in the  <dig> gene expression comparisons .
second, previous comparisons of bagged predictors in the context of genomic data were based on limited empirical evaluations.
many comparisons involved fewer than  <dig> microarray data sets when comparing predictors  <cit> .
while the comparisons involved clinically important data sets from cancer applications, these studies were simply not comprehensive enough.
third, previous studies probably did not consider enough bootstrap samples .
while previous studies used  <dig> to  <dig> bags, we always used  <dig> bags when constructing the rglm.
to illustrate how prediction accuracy depends on the number of bags, we evaluate the brain cancer data with  <dig> to  <dig> bags using  <dig> gene traits randomly selected from those used in our binary and continuous outcome prediction, respectively.
the results are shown in additional file  <dig>  most improvement is gained in the first several dozens of bags.
<dig> bags is generally enough although fluctuations remain.
more bags may lead to slightly better predictions but at the expense of longer computation time.
strengths and limitations
rglm shares many advantages of bagged predictors including a nearly unbiased estimate of the prediction accuracy  and several variable importance measures.
while our empirical studies focus on binary and continuous outcomes, it is straightforward to define rglm for count outcomes  and for multi-class outcomes .
a noteworthy limitation of rglm is computational complexity since the forward selection process  is particularly time-consuming.
the total time depends on the number of candidate features, the order of interaction terms, and the number of bags.
our r implementation allows the user to use parallel processing for speeding up the calculations.
our empirical studies demonstrate that rglm compares favorably with the random forest, support vector machines, penalized regression models, and many other widely used prediction methods.
as a caveat, we mention that we chose default parameter choices for each of these methods in order to ensure a fair comparison.
future studies could evaluate how these prediction methods compare when resampling schemes  are used to inform parameter choices.
our randomglm r package will allow the reader to carefully evaluate the method.
CONCLUSIONS
since individual forward selected glms are highly interpretable, the resulting ensemble predictor is more interpretable than an rf predictor.
our empirical studies  clearly highlight the outstanding prediction accuracy afforded by the rglm.
high accuracies are achieved not only in genomic data sets  but also in the uci benchmark data .
abbreviations
rglm: random generalized linear model; rglm: inter <dig> - rglm considering pairwise interactions between features; rglm: inter <dig> - rglm considering two-way and three-way interactions between features; forwardglm: forward selected generalized linear model; rf: random forest with default mtry; rfbigmtry: random forest with mtry equal to the total number of features; glm: generalized linear model; rpart: recursive partitioning; lda: linear discriminant analysis; dlda: diagonal linear discriminant analysis; knn: k nearest neighbor; svm: support vector machine; sc: shrunken centroids; rsm: random subspace method; rmnl: random multinomial logit model; rknn: random k nearest neighbor; e-rfe: entropy-based recursive feature elimination; aic: akaike information criteria; amv: adjusted majority vote.
competing interests
the authors declare that they have no competing interest.
authors’ contributions
ls carried out all analyses.
pl helped with the r implementation and analysis.
ls and sh developed the method and wrote the article.
sh conceived of the study.
all authors read and approved the final manuscript.
supplementary material
additional file 1
simulation study design.
this file describes the simulation studies and presents r code used for simulating the data set.
click here for file
additional file 2
sensitivity and specificity of predictors in the  <dig> disease gene expression data sets.
for each data set and prediction method, the table reports the sensitivity and specificity estimated using 3-fold cross validation.
more precisely, the table reports the average 3-fold cv estimate over  <dig> random partitions of the data into  <dig> folds.
median sensitivity and specificity across data sets are summarized at the bottom.
click here for file
additional file 3
sensitivity and specificity of predictors in the uci machine learning benchmark data.
for each data set and prediction method, the table reports the sensitivity and specificity estimated using 3-fold cross validation.
more precisely, the table reports the average 3-fold cv estimate over  <dig> random partitions of the data into  <dig> folds.
median sensitivity and specificity across data sets are summarized at the bottom.
click here for file
additional file 4
prediction accuracy when including pairwise interactions between features in the uci machine learning benchmark data.
this table is an extension to table  <dig>  it shows the prediction accuracy of predictors other than rglm when considering pairwise interactions between features in the same uci mlbench data sets.
although several predictors show improvement, none of them beats rglm.inter <dig>
click here for file
additional file 5
comparison of rglm based feature selection method with the rf based method of díaz-uriarte et al.
for each data set in the  <dig> disease gene expression data, the rf based variable selection method by díaz-uriarte et al selects a small set of genes.
for each of the selected genes, the file reports the ranking in terms of the rglm variable importance measure timesselectedbyforwardregression.
as expected, only a few of the selected genes have a high rank in terms of timesselectedbyforwardregression illustrating that these variable selection methods are different.
click here for file
additional file 6
effect of the number of bags on rglm predictor thinning.
s this figure reports how prediction accuracy changes as variable thinning is applied to the rglm.
results are averaged over the  <dig> dichotomized gene traits in the mouse adipose data set.
the five rows correspond to nbags values of  <dig>   <dig>   <dig>   <dig>   <dig> respectively.
within each row, the two panels have the same meaning as in figure  <dig>
click here for file
additional file 7
prediction accuracy versus number of bags used for rglm.
this figure presents the results for predicting  <dig> gene traits in the brain cancer data set when different numbers of bags  are used for constructing the rglm.
each color represents one gene trait.
binary outcome prediction.
the  <dig> gene traits were randomly selected from all  <dig> gene traits used in the binary outcome prediction section.
continuous outcome prediction.
the  <dig> gene traits were randomly selected from all  <dig> gene traits used in the continuous outcome prediction.
click here for file
acknowledgements
we acknowledge grant support from 1r01da030913- <dig>  p50ca <dig>  p30ca <dig>  ul1tr <dig>  we acknowledge the efforts of igc and expo in providing data set gse <dig>

 prediction methods are widely used in biomedical research. for example, reliable prediction methods are essential for accurate disease classification, diagnosis and prognosis. since prediction methods based on multiple features can greatly outperform predictors based on a single feature, it is important to develop methods that can optimally combine features to obtain high accuracy. introductory text books describe well known prediction methods such as linear discriminant analysis, k-nearest neighbor predictors, support vector machines, and tree predictors. many publications have evaluated popular prediction methods in the context of gene expression data. ensemble predictors are particularly attractive since they are known to lead to highly accurate predictions. an ensemble predictor generates and integrates multiple versions of a single predictor, and arrives at a final prediction by aggregating the predictions of multiple base learners, e.g. via plurality voting across the ensemble. one particular approach for constructing an ensemble predictor is bootstrap aggregation. here multiple versions of the original data are generated through bootstrapping, where observations from the training set are randomly sampled with replacement. an individual predictor is fitted on each bootstrapped data set. thus, bootstrapped data sets will lead to an ensemble of tree predictors. in case of a class outcome, the individual predictors “vote” for each class and the final prediction is obtained by majority voting. breiman showed that bagging weak predictors often yields substantial gains in predictive accuracy. but it seems that ensemble predictors are only very rarely used for predicting clinical outcomes. this fact points to a major weakness of ensemble predictors: they typically lead to ”black box” predictions that are hard to interpret in terms of the underlying features. clinicians and epidemiologists prefer forward selected regression models since the resulting predictors are highly interpretable: a linear combination of relatively few features can be used to predict the outcome or the probability of an outcome. but the sparsity afforded by forward feature selection comes at an unacceptably high cost: forward variable selection often greatly overfit the data which results in unstable and inaccurate predictors. ideally, one would want to combine the advantages of ensemble predictors with those of forward selected regression models. as discussed below, multiple articles describe ensemble predictors based on linear models including the seminal work by breiman who evaluated a bagged forward selected linear regression model. however, the idea of bagging forward selected linear models appears to have been set aside as new ensemble predictors, such as the random forest, became popular. a random forest predictor not only bags tree predictors but also introduces an element of randomness by considering only a randomly selected subset of features at each node split. the number of randomly selected features, mtry, is the only parameter of the random forest predictor. the random forest predictor has deservedly received a lot of attention for the following reasons: first, the bootstrap aggregation step allows one to use out-of-bag samples to estimate the predictive accuracy. the resulting oob estimate of the accuracy often obviates the need for cross-validation and other resampling techniques. second, the rf predictor provides several measures of feature importance. several articles explore the use of these importance measures to select genes. third, it can be used to define a dissimilarity measure that can be used in clustering applications. fourth, and most importantly, the rf predictor has superior predictive accuracy. it performs as well as alternatives in cancer gene expression data sets but it really stands out when applied to the uci machine learning benchmark data sets where it is as good as many existing methods. while we confirm the truly outstanding predictive performance of the rf, the proposed rglm method turns out to be even more accurate than the rf. breiman and others have pointed out that the black box predictions of the rf predictor can be difficult to interpret. for this reason, we wanted to give bagged forward selected generalized linear regression models another careful look. after exploring different approaches for injecting elements of randomness into the individual glm predictors, we arrived at a new ensemble predictor, referred to as random glm predictor, with an astonishing predictive performance. an attractive aspect of the proposed rglm predictor is that it combines the advantages of the rf with that of a forward selected glm. as the name generalized linear model indicates, it can be used for a general outcome such as a binary outcome, a multi-class outcome, a count outcome, and a quantitative outcome. we show that several incremental changes to the original bagged glm predictor by breiman add to up to a qualitatively new predictor that performs at least as well as the rf predictor on the uci benchmark data sets. while the uci data are the benchmark data for evaluating predictors, only a dozen such data sets are available for binary outcome prediction. to provide a more comprehensive empirical comparison of the different prediction methods, we also consider over comparisons involving gene expression data. in these genomic studies, the rglm method turns out to be slightly more accurate than the considered alternatives. while the improvements in accuracy afforded by the rglm are relatively small they are statistically significant. this article is organized as follows. first, we present a motivating example that illustrates the high prediction accuracy of the rglm. second, we compare the rglm with other state of the art predictors when it comes to binary outcome prediction. toward this end, we use the uci machine learning benchmark data, over empirical gene expression comparisons, and extensive simulations. third, we compare the rglm with other predictors for quantitative outcome prediction. fourth, we describe several variable importance measures and show how they can be used to define a thinned version of the rglm that only uses few important features. even for data sets comprised of thousands of gene features, the thinned rglm often involves fewer than features and is thus more interpretable than most ensemble predictors. methods construction of the rglm predictor rglm is an ensemble predictor based on bootstrap aggregation of generalized linear models whose features are selected using forward regression according to aic criterion. glms comprise a large class of regression models, e.g. linear regression for a normally distributed outcome, logistic regression for binary outcome, multi-nomial regression for multi-class outcome and poisson regression for count outcome. thus, rglm can be used to predict binary-, continuous-, count-, and other outcomes for which generalized linear models can be defined. the “randomness” in rglm stems results from two sources. first, a non-parametric bootstrap procedure is used which randomly selects samples with replacement from the original data set. second, a random subset of features is selected for each bootstrap sample. this amounts to a random sub-space method applied to each bootstrap sample separately. the steps of the rglm construction are presented in figure first, starting from the original data set another equal-sized data set is generated using the non-parametric bootstrap method, i.e. samples are selected with replacement from the original data set. the parameter nbags determines how many of such bootstrap data sets are being generated. second, a random set of features is randomly chosen for each bag. thus, the glm predictor for bag will typically involve a different set of features than that for bag third, the nfeaturesinbag of randomly selected features per bag are rank-ordered according to their individual association with the outcome variable y in each bag. for a quantitative outcome y, one can simply use the absolute value of the correlation coefficient between the outcome and each feature to rank the features. more generally, one can fit a univariate glm model to each feature to arrive at an association measure. only the top ranking features will become candidate covariates for forward selection in a multivariate regression model. the top number of candidate features is determined by the input parameter ncandidatecovariates. fourth, forward variable selection is applied to the ncandidatecovariates of each bag to arrive at a multivariable generalized linear model per bag. the forward selection procedure used by rglm is based on the stepaic r function in the mass r library where method is set to “forward”. fifth, the predictions of each forward selected multivariate model are aggregated across bags to arrive at a final ensemble prediction. the aggregation method depends on the type of outcome. for a continuous outcome, predicted values are simply averaged across bags. for a binary outcome, the adjusted majority vote strategy is used which averages predicted probabilities across bags. given the estimated class probabilities one can get a binary prediction by choosing an appropriate threshold. importantly, rglm also has a parameter maxinteractionorder for creating interactions up to a given order among features in the model construction. for example, rglm.inter results from setting maxinteractionorder= i.e. considering pairwise interaction terms. as example, consider the case when only pairwise interaction terms are used. for each bag a random set of features is selected from the original covariates, i.e. covariates without interaction terms. next, all pairwise interactions among the nfeaturesinbag randomly selected features are generated. next, the usual rglm candidate feature selection steps will be applied to the combined set of pairwise interaction terms and the nfeaturesinbag randomly selected features per bag resulting in ncandidatecovariates top ranking features per bag, which are subsequently subjected to forward feature selection. these methods are implemented in our r software package randomglm which allows the user to input a training set and optionally a test set. it automatically outputs out-of-bag estimates of the accuracy and variable importance measures. parameter choices for the rglm predictor as discussed below, we find that it is usually sufficient to consider only nbags= bootstrap data sets. the default value for nfeaturesinbag depends on the total number of features. it is easier to explain it in terms of the proportion of features randomly selected per bag, nfeaturesinbag/n, where n is the total number of features of the training set. apart from n it is also valuable to consider the effective number of features which equals the number of features n plus the number of interaction terms, e.g. n∗=n+n/ in case of pairwise interactions. using this notation, the default value of nfeaturesinbag can be arrived at by solving equations presented in table these equations were found by empirically evaluating various choices of nfeaturesinbag values. in particular, we found that in case of n∗<= then using all features is often a good choice, whereas if n∗> then setting nfeaturesinbag/n= works well. the default value nfeaturesinbag/n= − n∗ in the intermediate case results from fitting an interpolation line through the two points and. we find that rglm is quite robust with respect to the parameter nfeaturesinbag. to limit the number of covariates considered in forward selection, the default value of ncandidatecovariates is set to overall, the default values perform well in our simulations, empirical gene expression and machine learning benchmark studies. but we recommend to use the oob estimate of predictive accuracy to inform the choice of the parameter values. this table shows the default values of nfeaturesinbag in terms of nfeaturesinbag/n for rglm, rglm.inter and rglm.inter n is the total number of features of the training data. n∗ is the effective number of features which equals the number of features n plus the number of interaction terms. formulas are shown in terms of both n and the corresponding n∗. and are obtained by interpolating a straight line between and. relationship with related prediction methods as discussed below, rglm can be interpreted as a variant of a bagged predictor. in particular, it is similar to the bagged forward linear regression model but differs in the following aspects: rglm allows for interaction terms between features which greatly improve the performance on some data sets. we refer to rglm involving two-way or three way interactions as rglm.inter and rglm.inter respectively. rglm has a parameter nfeaturesinbag that allows one to restrict the number of features used in each bootstrap sample. this parameter is conceptually related to the mtry parameter of the random forest predictor. in essence, this parameter allows one to use a random subspace method in each bootstrap sample. rglm has a parameter ncandidatecovariates that allows one to restrict the number of features in forward regression, which not only has computational advantages but also introduces additional instability into the individual predictors, which is a desirable characteristic of an ensemble predictor. rglm optimizes the aic criterion during forward selection. rglm has a “thinning threshold” parameter which allows one to reduce the number of features involved in prediction while maintaining good prediction accuracy. since a thinned rglm involves far fewer features, it facilitates the understanding how the ensemble arrives at its predictions. rglm is not only related to bagging but also to the random subspace method proposed by. in the rsm, the training set is also repeatedly modified as in bagging but this modification is performed in the feature space. in the rsm, a subset of features is randomly selected which amounts to restricting attention to a subspace of the original feature space. as one of its construction steps, rglm uses a rsm on each bootstrap sample. future research could explore whether random partitions as opposed to random subspaces would be useful for constructing an rglm. random partitions of the feature space are similar to random subspaces but they divide the feature space into mutually exclusive subspaces. random partition based predictors have been shown to perform well in high-dimensional data. both rsm and random partitions have more general applicability than rglm since these methods can be used for any base learner. there is a vast literature on ensemble induction methods but a property worth highlighting is that rglm uses forward variable selection of glms. recall that rglm goes through the following steps: 1) bootstrap sampling, 2) rsm, 3) forward variable selection of a glm, 4) aggregation of votes. empirical studies involving different base learners have shown that combining bootstrap sampling with rsm leads to ensemble predictors with comparable performance to that of the random forest predictor. another prediction method, random multinomial logit model, also shares a similar idea with rglm. it was recently proposed for multi-class outcome prediction. rmnl bags multinomial logit models with random feature selection in each bag. it can be seen as a special case of rglm, except that no forward model selection is carried out. software implementation the rglm method is implemented in the freely available r package randomglm. the r function randomglm allows the user to output training set predictions, out-of-bag predictions, test set predictions, coefficient values, and variable importance measures. the predict function can be used arrive at test set predictions. tutorials can be found at the following webpage: http://labs.genetics.ucla.edu/horvath/rglm. short description of alternative prediction methods forward selected generalized linear model predictor we denote by forwardglm the generalized linear model predictor whose covariates were selected using forward feature selection. thus, forwardglm does not involve bagging, random feature selection, and is not an ensemble predictor. random forest rf is an ensemble predictor that consists of a collection of decision trees which vote for the class of observations. the rf is known for its outstanding predictive accuracy. we used the randomforest r package in our studies. we considered two choices for the rf parameter mtry: i) the default rf predictor where mtry equals the square root of the number of features and ii) rfbigmtry where mtry equals the total number of features. we always generated at least trees per forest but used trees when calculating variable importance measures. recursive partitioning and regression trees classification and regression trees were generated using the default settings rpart r package. tree methods are described in. linear discriminant analysis lda aims to find a linear combination of features to predict a binary outcome. we used the lda r function in the mass r package with parameter choice method=moment. diagonal linear discriminant analysis dlda is similar to lda but it ignores the correlation patterns between features. while this is often an unrealistic assumption, dlda has been found to work well in in gene expression applications. here we used the default parameters from the supclust r package. k nearest neighbor we used the knn r function in the class r package, which chose the parameter k of nearest neighbors using 3-fold cross validation. support vector machines we used the default parameters from the e r package to fit svms. additional details can be found in. shrunken centroids the sc predictor is known to work well in the context of gene expression data. here we used the implementation in the pamr r package which chose the optimal level of shrinkage using cross validation. penalized regression models various convex penalties can be applied to generalized linear models. we considered ridge regression corresponding to an ℓ penalty, the lasso corresponding to an ℓ penalty, and elastic net corresponding to a linear combination of ℓ and ℓ penalties. we used the glmnet r function from the glmnet r package with alpha parameter values of and respectively. glmnet also involves another parameter which was chosen as the median of the lambda sequence output resulting from glmnet. for uci benchmark data sets, pairwise interaction between features were considered. disease-related gene expression data sets we use disease related gene expression data sets involving cancer and other human diseases. the first data sets involving various cancers were previously used by. these data can be downloaded from the author’s webpage at http://ligarto.org/rdiaz/papers/rfvs/randomforestvarsel.html. the braintumor and dlbcl data sets were downloaded from http://www.gems-system.org/. the remaining data sets were downloaded from either the gene expression omnibus database or the arrayexpress data base in raw form and subsequently preprocessed using mas normalization and quantile normalization. only the top probes with highest mean expression were considered for outcome prediction. we briefly point out that diaz et al report prediction error rates estimated using a bootstrap method. in contrast, we report 3-fold cross validation estimates, which may explain minor numerical differences between our study and that of diaz et al. sample size, number of features, original reference, data set ids and outcomes for the disease related gene expression data sets. empirical gene expression data sets for all data sets below, we considered randomly selected gene traits, i.e. randomly selected probes. they were directly used as continuous outcomes or dichotomized according to the median value to generate binary outcomes. for all data sets except “brain cancer”, of the observations were randomly chosen as the training set, while the remaining samples were chosen as test set. we focused on the genes with the highest mean expression levels in each data set. brain cancer data sets these two related data sets contain and microarray samples of glioblastoma patients, respectively. gene expression profiles were measured using affymetrix u microarrays. a detailed description can be found in. the first data set was used as a training set while and the second data set was used as a test set. safhs blood lymphocyte data set this data set was derived from blood lymphocytes of randomly ascertained participants enrolled in the san antonio family heart study. gene expression profiles were measured with the illumina sentrix human whole genome series i beadchips. after removing potential outliers, samples remained in the data set. wb whole blood gene expression data set this is the whole blood gene expression data from healthy controls. peripheral blood samples from healthy individuals were analyzed using illumina human ht- microarrays. after pre-processing, samples remained in the data set. mouse tissue gene expression data sets the tissue specific gene expression data sets were generated by the lab of jake lusis at ucla. these data sets measure gene expression levels from adipose, brain, liver and muscle tissue of mice from the b ×h f mouse intercross described in. in addition to gene traits, we also predicted quantitative mouse clinical traits including mouse weight, length, abdominal fat, other fat, total fat, adiposity index, plasma triglycerides, total plasma cholesterol, high-density lipoprotein fraction of cholesterol, plasma unesterified cholesterol, plasma free fatty acids, plasma glucose, plasma low-density lipoprotein and very low-density lipoprotein cholesterol, plasma mcp- protein levels, plasma insulin, plasma glucose-insulin ratio, plasma leptin, plasma adiponectin, aortic lesion size, aneurysms, and aortic calcification in the lesion area. machine learning benchmark data sets the machine learning benchmark data sets used in this article are listed in table note that only eight of the data sets have a binary outcomes. the multi-class outcomes of the remaining data sets were turned into binary outcomes by considering the most prevalent class versus all other classes combined. missing data were imputed using nearest neighbor averaging. for each data set and prediction method, we report the average 3-fold cv estimate of prediction accuracy over random partitions of the data into folds. sample size and number of features for the uci machine learning benchmark data sets. simulated gene expression data sets we simulated an outcome variable y and gene expression data that contained modules. only of the modules were comprised of genes that correlated with the outcome y. 45% of the genes were background genes, i.e. these genes were outside of any module. the simulation scheme is detailed in additional file and implemented in the r function simulatedatexpr5modules from the wgcna r package. this r function was used to simulate pairs of training and test data sets. the simulation study was used to evaluate prediction methods for continuous outcomes and for binary outcomes. for binary outcome prediction, the continuous outcome y was thresholded according to its median value. we considered different simulation scenarios involving varying sizes of the training data and varying numbers of genes that served as features. test sets contained the same number of genes as in the corresponding training set and samples. for each simulation scenario, we simulate replicates resulting from different choices of the random seed.  motivating example: disease-related gene expression data sets we compare the prediction accuracy of rglm with that of other widely used methods on gene expression data sets involving human disease related outcomes. many of the data sets are well known cancer data sets, which have been used in other comparative studies. a brief description of the data sets can be found in methods. to arrive at an unbiased estimate of prediction accuracy, we used 3-fold cross validation. note that the accuracy equals minus the median misclassification error rate. table reports the prediction accuracy of different methods including rglm, random forest, random forest, tree predictor, linear discriminant analysis, diagonal linear discriminant analysis, k nearest neighbor, support vector machine and shrunken centroid. a short description of these prediction methods is provided in methods. for each data set, the prediction accuracy was estimated using 3−fold cross validation across random partitions of the data into folds. mean accuracies across the data sets and the resulting ranks are summarized at the bottom. two sided paired wilcoxon test p-values can be used to determine whether the accuracy of rglm is significantly different from that of other predictors. note that the rglm yields the highest mean accuracy. as seen from table rglm achieves the highest mean accuracy in these disease data sets, followed by rfbigmtry and sc. note that the standard random forest predictor performs worse than rglm. the accuracy difference between rglm and alternative methods is statistically significant for all predictors except for rfbigmtry, dlda and sc. since rfbigmtry is an ensemble predictor that relies on thousands of features it would be difficult to interpret its predictions in terms of the underlying genes. our evaluations focused on the accuracy. however, a host of other accuracy measures could be considered. additional file presents the results for sensitivity and specificity. the top methods with highest sensitivity are: rf, svm and rglm. the top methods with highest specificity are: sc, rglm and knn. a strength of this empirical comparison is that it involves clinically or biologically interesting data sets but a severe limitation is that it only involves comparisons. therefore, we now turn to more comprehensive empirical comparisons. binary outcome prediction empirical study involving dichotomized gene traits many previous empirical comparisons of gene expression data considered fewer than data sets. to arrive at comparisons, we use the following approach: we started out with human and mouse gene expression data sets. for each data set, we randomly chose genes as gene traits resulting in 7× possible outcomes. we removed the gene corresponding to the gene trait from the feature set. next, each gene trait was dichotomized by its median value to arrive at a binary outcome y. the goal of each prediction analysis was to predict the dichotomized gene trait y based on the other genes. at first sight, this artificial outcome is clinically uninteresting but it is worth emphasizing that clinicians often deal with dichotomized measures of gene products, e.g. high serum creatinine levels may indicate kidney damage, high psa levels may indicate prostate cancer, and high hdl levels may indicate hypercholesterolemia. to arrive at unbiased estimates of prediction accuracy, we split each data set into a training and test set. figure shows boxplots of the accuracies across the comparisons. similar performance patterns are observed for the individual data sets ). the figure also reports pairwise comparisons of the rglm method versus alternative methods. specifically, it reports the two-sided wilcoxon signed rank test p-values for testing whether the accuracy of the rglm predictor is higher than that of the considered alternative method. strikingly, rglm is more accurate than the other methods overall. while the increase in accuracy is often minor, it is statistically significant as can be seen by comparing rglm to rf, rfbigmtry, lda, svm and sc. other predictors perform even worse, and the corresponding p-values are not shown. the fact that rfbigmtry is more accurate in this situation than the default version of rf probably indicates that relatively few genes are informative for predicting a dichotomized gene trait. also note that rglm is much more accurate than the unbagged forward selected glm which reflects that forward selection greatly overfits the training data. in conclusion, these comprehensive gene expression studies show that rglm has outstanding prediction accuracy. machine learning benchmark data analysis here we evaluate the performance of rglm on the uci machine learning benchmark data sets which are often used for evaluating prediction methods. we consider benchmark data sets from the mlbench r package: uci data sets and synthetic data sets. we choose these data sets for two reasons. first, these data sets were also used in the original evaluation of the random forest predictor. second, these data include all of the available data sets with binary outcomes in the mlbench r package. a detailed description of these data sets can be found in methods. in his original publication on the random forest, breiman found that the rf outperformed bagged predictors on the uci benchmark data which may explain why bagged glms have not received much attention. we hypothesize that the relatively poor performance of a bagged logistic regression model on these data sets could be ameliorated by considering interaction terms between the features. table confirms our hypothesis. rglm.inter has superior or tied accuracy compared to rglm in out of benchmark data sets. in particular, pairwise interactions greatly improve the prediction accuracy in the ringnorm data set. higher order interactions do not perform better than rglm.inter but dramatically increase computational burden. for each data set, the prediction accuracy was estimated using 3−fold cross validation across random partitions of the data into folds. rglm.inter incorporates pairwise interaction between features into the rglm predictor. mean accuracies and the resulting ranks are summarized at the bottom. the wilcoxon signed rank test was used to test whether accuracy differences between rglm.inter and other predictors are significant. rglm.inter rf, and svm tie for first place. overall, we find that rglm.inter ties with svm and rf for the first place in the benchmark data. as can be seen from additional file rglm.inter achieves the highest sensitivity and specificity, which also support its good performance in the benchmark data sets. a potential limitation of these comparisons is that we considered pairwise interaction terms for the rglm predictor but not for the other predictors. to address this issue, we also considered pairwise interactions among features for other predictors. additional file shows that no method surpasses rglm.inter when pairwise interaction terms are considered. in particular, interaction terms between features do not improve the performance of the random forest predictor. a noteworthy disadvantage of rglm.inter in case of many features is the computational burden that may result from adding interaction terms. in applications where interaction terms are needed for rglm, faster alternatives remain an attractive choice. simulation study involving binary outcomes as described in methods, we simulated gene expression data sets with binary outcomes. the number of features ranged from to the sample sizes of the training data ranged from to to robustly estimate the test set accuracy we chose a large size for the corresponding test set data, n= figure shows the boxplots of the test set accuracies of different predictors. the accuracy of the forwardglm is much lower than that of rglm, demonstrating the benefit of creating an ensemble predictor. the wilcoxon test p-value shows that rglm is significantly better than all other methods except for the rf. in this simulation study, rglm takes the first place when it comes to the median accuracy. continuous outcome prediction in the following, we show that rglm also performs exceptionally well when dealing with continuous quantitative outcomes. we not only compare rglm to a standard forward selected linear model predictor but also a random forest predictor. we do not report the findings for the k-nearest neighbor predictor of a continuous outcome since it performed much worse than the above mentioned approaches in our gene expression applications. we again split the data into training and test sets. we use the correlation between test set predictions and truly observed test set outcomes as measure of predictive accuracy. note that this correlation coefficient can take on negative values. empirical study involving continuous gene traits here we used the same gene expression comparisons as described above but did not dichotomize the gene traits. incidentally, prediction methods for gene traits are often used for imputing missing gene expression values. our results presented in figure indicate that for the majority of genes high accuracies can be achieved. but for some gene traits, the accuracy measure, which is defined as a correlation coefficient, takes on negative values indicating that there is no signal in the data. note that the forward selected linear predictor ties with the random forest irrespective of the choice of the mtry parameter and both methods perform significantly worse than the rglm predictor. mouse tissue expression data involving continuous clinical outcomes here we used the mouse liver and adipose tissue gene expression data sets to predict clinical outcomes. again, rglm achieved significantly higher median prediction accuracy compared to the other predictors. simulation study involving continuous outcomes gene expression data sets are simulated in the same way as described previously but here the outcome y was not dichotomized. as shown in figure the forwardglm accuracy trails both rglm and rf, reflecting again the fact that forward regression overfits the data. in this simulation study, we find that rglm yields significantly higher prediction accuracy than other predictors. comparing rglm with penalized regression models in our previous comparisons, we found that rglm greatly outperforms forward selected glm methods based on the aic criterion. many powerful alternatives to forward variable selection have been developed in the literature, in particular penalized regression models. here, we compare rglm to major types of penalized regression models: ridge regression, elastic net, and the lasso. the predictive accuracies of these penalized regression models were compared to those of the rglm predictor using the same data sets described above for evaluating binary outcome and quantitative outcome prediction methods. wilcoxon’s signed rank test was used to determine whether differences in predictive accuracy were significant. figure shows that rglm outperforms penalized regression models when applied to binary outcomes. for all comparisons, the paired median difference is positive which indicates that rglm is at least as good if not better than any of these penalized regression models. in particular, rglm is significantly better than ridge regression and the lasso on the dichotomized gene expression trait data. also, rglm is significantly better than elastic net and lasso in simulations with binary outcomes. figure shows that rglm outperforms penalized regression models for continuous outcome prediction as well. positive accuracy differences again imply that rglm is at least as good as these penalized regression models. in particular, it significantly outperforms ridge regression in the continuous gene expression traits data and outperforms elastic net and lasso in simulations with continuous outcomes. as a caveat, we mention that cross validation methods were not used to inform the parameter choices of the penalized regression models since the rglm predictor was also not allowed to fine tune its parameters. by only using default parameter choices we ensure a fair comparison. in a secondary analysis, however, we allowed penalized regression models to use cross validation for informing the choice of the parameters. while this slightly improved the performance of the penalized regression models, it did not affect our main conclusion. rglm outperforms penalized regression models in these comparisons. feature selection here we briefly describe how rglm naturally gives rise to variable importance measures. we compare the variable importance measures of rglm with alternative approaches and show how variable importance measures can be used for defining a thinned rglm predictor with few features. variable importance measure there is a vast literature on using ensemble predictors and bagging for selecting features. for example, meinshausen and bühlmann describe “stability selection” based on variable selection employed in regression models. the method involves repetitive sub-sampling, and variables that occur in a large fraction of the resulting selection set are chosen. li et al. use a random k-nearest neighbor predictor to carry out feature selection. the entropy-based recursive feature elimination method of furlanello et al. ranks features in high dimensional microarray data. rglm, like many ensemble predictors, gives rise to several measures of feature importance. for example, the number of times a feature is selected in the forward glm across bags, timesselectedbyforwardregression, is a natural measure of variable importance. another variable importance measure is the number of times a feature is selected as candidate covariate for forward regression, timesselectedascandidates. note that both timesselectedbyforwardregression and timesselectedascandidates have to be ≤nbags. finally, one can use the sum of absolute glm coefficient values, sumabscoefbyforwardregression, as a variable importance measure. we prefer timesselectedbyforwardregression, since it is more intuitive and points to the features that directly contribute to outcome prediction. to reveal relationships between different types of variable importance measures, we present a hierarchical cluster tree of rglm measures, rf measures and standard marginal analysis based on correlations in figure as expected, the marginal association measures cluster together. the same holds for the random forest based importance measures and the rglm based importance measures. leo breiman already pointed out that random forests could be used for feature selection in genomic applications. díaz-uriarte et al. proposed a related gene selection method based on the rf which yields small sets of genes. this rf based gene selection method does not return sets of genes that are highly correlated because such genes would be redundant when it comes to predictive purposes. since the rglm based importance measure timesselectedbyforwardregression is expected to lean towards selecting genes that are highly associated with the outcome, it comes as no surprise that only a few genes selected by the procedure of díaz-uriarte et al. turn out to have a top ranking in terms of the rglm measure timesselectedbyforwardregression. it is beyond our scope to provide a thorough evaluation of the different variable selection approaches and we refer the reader to the literature, e.g.. while our studies show that the rglm based variable importance measures have some relationships to other measures, they are sufficiently different from other measures to warrant a thorough evaluation in future comparison studies. rglm predictor thinning based on a variable importance measure both rglm and random forest have superior prediction accuracy but they differ with respect to how many features are being used. recall that the random forest is composed of individual trees. each tree is constructed by repeated node splits. the number of features considered at each node split is determined by the rf parameter mtry. the default value of mtry is the square root of the number of features. in case of gene features in our empirical studies, the default value is mtry= for rfbigmtry, we choose all possible features, i.e. mtry= we find that a random forest predictor typically uses more than 40% of the features in the empirical studies. in contrast, rglm typically only involves a few hundred genes in these studies. there are several reasons why rglm uses far fewer features in its construction. first, and foremost, it uses forward selection to select features in each bag. second, the number of candidate covariates considered for forward regression is chosen to be low, i.e. ncandidatecovariates= in rglm, the number of times a feature is selected by forward regression models among all bags, timesselectedbyforwardregression, follows a highly skewed distribution. only few features are repeatedly selected into the model while most features are selected only once. it stands to reason that an even sparser, highly accurate predictor can be defined by refitting the glm on each bag without considering these rarely selected features. we refer to this feature removal process as rglm predictor thinning. thus, features whose value of timesselectedbyforwardregression lies below a pre-specified thinning threshold will be removed from the model fit a posteriori. interestingly, the accuracy diminishes very slowly for initial, low threshold values. but even low threshold values lead to a markedly sparser ensemble predictor ). in other words, the average fraction of features remaining in the thinned rglm declines drastically as the thinning threshold increases. we have found that the following empirical function accurately describes the relationship between thinning threshold and proportion of features left in the thinned rglm predictor: proportionleft=f=1x=0exp{−e nbags )}0<x≤ where x=thinningthresholdnbags and e denotes euler’s constant e≈. equation was found by log transforming the data and using optimization approaches for estimating the parameters. no mathematical derivation was used. one can easily show that f is a monotonically decreasing function which accurately describes the proportion of remaining features as can be seen from figure. since the proportion of remaining variables depends not only on the thinning threshold but also on the number of bags nbags, we also study how these results depend on the choice of nbags. toward this end, we varied nbags from to for predicting the dichotomized gene traits in the mouse adipose data set. additional file shows that the predicted values based on equation overlaps almost perfectly with the observed values for all considered choices of nbags, which indicates that equation accurately estimates the proportion of remaining features for range of different values of nbags. our results demonstrate that the number of required features decreases rapidly even for low values of the thinning threshold without compromising the prediction accuracy of the thinned predictor. figure shows that a thinning threshold of leads to a thinned predictor whose accuracy is negligibly lower than that of the original rglm predictor but it involves less than 20% of the original number of variables. recall that even the original number of variables is markedly lower than that of the rf predictor. these results demonstrate that the thinned rglm combines the advantages of an ensemble predictor with that of a forward selected glm model. rglm thinning versus rf thinning the idea behind rglm thinning is to remove features with low values of the variable importance measure. of course, a similar idea can be applied to other predictors. here we briefly evaluate the performance of a thinned random forest predictor which removed variables based on a low value of its importance measure. to arrive at an unbiased comparison, both rglm and rf are thinned based on results obtained in the training data. next, accuracies of the thinned predictors are evaluated in the test set data. figure compares thinned rglm versus thinned rf in our disease related data sets and also the empirical studies. numbers that connect dashed lines are rglm thinning thresholds. for a pre-specified threshold, the number of features used in the thinned random forest is matched to that used in the thinned rglm. without thinning, rf uses a lot more features than rglm as mentioned previously. as expected, the median number of genes left for prediction and the corresponding median prediction accuracy generally decrease as the thinning threshold becomes more stringent. overall, a thinned rglm yields a significantly higher median accuracy than a thinned rf across different thinning thresholds. in clinical practice, a thinned predictor with very few features and good accuracy can be very useful and interpretable. for example, choosing a threshold of in panel figure and a threshold of in panel would result in very sparse predictors. in both cases, especially in panel, the thinned rglm has higher median accuracy than that of the thinned rf. discussion why was the rglm not discovered earlier? after breiman proposed the idea of bagged linear regression models in, many authors have explored the utility of bagging logistic regression models. most previous studies report that bagging does not improve the accuracy of logistic regression. bühlman and yu showed theoretically that bagging helps for “hard threshold” methods but not for “soft threshold” methods. these studies indicate that bagged logistic regression models are not beneficial since the individual predictors are too stable. overall, we agree with these results. but our comprehensive evaluations show that by injecting elements of randomness and instability into a bagged logistic regression model one arrives at a state of the art prediction method that often outperforms existing methods. figure describes why the construction of the rglm runs counter to conventional wisdom. as indicated by the upper right hand panel of figure the rglm is based on two seemingly bad modifications to a glm. as indicated by the top left panel of figure forward selection of a glm is typically a bad idea since it overfits the data and thus degrades the prediction accuracy of a single glm predictor. as indicated by the bottom right panel of figure bagging a full logistic regression is also a bad idea since it leads to a complicated predictor without clear evidence for increased accuracy. but these two seemingly bad modifications add up to a superior prediction method. breiman already noted that the instability afforded by variable selection is important for constructing a bagged linear model based predictor. in order to define an accurate glm based ensemble predictor, we also find that it is important to introduce additional elements of randomness and instability, which is also reflected in the name random glm. our results show that the proposed changes results in a more accurate predictor that involves surprisingly few features. additional reasons why the merits of rglm have not been recognized earlier may be the following. first, it may be a historical accident. bagging was quickly over-shadowed by other seemingly more accurate ways of constructing ensemble predictors, such as boosting and the rf, both of which have markedly better performance on the uci benchmark data. we find that rglm.inter ties with svm and rf for the top spot in uci benchmark data set. incidentally, rglm performs significantly better than svm and rf on the disease data sets and in the gene expression comparisons. second, previous comparisons of bagged predictors in the context of genomic data were based on limited empirical evaluations. many comparisons involved fewer than microarray data sets when comparing predictors. while the comparisons involved clinically important data sets from cancer applications, these studies were simply not comprehensive enough. third, previous studies probably did not consider enough bootstrap samples. while previous studies used to bags, we always used bags when constructing the rglm. to illustrate how prediction accuracy depends on the number of bags, we evaluate the brain cancer data with to bags using gene traits randomly selected from those used in our binary and continuous outcome prediction, respectively. the results are shown in additional file most improvement is gained in the first several dozens of bags. bags is generally enough although fluctuations remain. more bags may lead to slightly better predictions but at the expense of longer computation time. strengths and limitations rglm shares many advantages of bagged predictors including a nearly unbiased estimate of the prediction accuracy and several variable importance measures. while our empirical studies focus on binary and continuous outcomes, it is straightforward to define rglm for count outcomes and for multi-class outcomes. a noteworthy limitation of rglm is computational complexity since the forward selection process is particularly time-consuming. the total time depends on the number of candidate features, the order of interaction terms, and the number of bags. our r implementation allows the user to use parallel processing for speeding up the calculations. our empirical studies demonstrate that rglm compares favorably with the random forest, support vector machines, penalized regression models, and many other widely used prediction methods. as a caveat, we mention that we chose default parameter choices for each of these methods in order to ensure a fair comparison. future studies could evaluate how these prediction methods compare when resampling schemes are used to inform parameter choices. our randomglm r package will allow the reader to carefully evaluate the method.  since individual forward selected glms are highly interpretable, the resulting ensemble predictor is more interpretable than an rf predictor. our empirical studies clearly highlight the outstanding prediction accuracy afforded by the rglm. high accuracies are achieved not only in genomic data sets but also in the uci benchmark data. abbreviations rglm: random generalized linear model; rglm: inter - rglm considering pairwise interactions between features; rglm: inter - rglm considering two-way and three-way interactions between features; forwardglm: forward selected generalized linear model; rf: random forest with default mtry; rfbigmtry: random forest with mtry equal to the total number of features; glm: generalized linear model; rpart: recursive partitioning; lda: linear discriminant analysis; dlda: diagonal linear discriminant analysis; knn: k nearest neighbor; svm: support vector machine; sc: shrunken centroids; rsm: random subspace method; rmnl: random multinomial logit model; rknn: random k nearest neighbor; e-rfe: entropy-based recursive feature elimination; aic: akaike information criteria; amv: adjusted majority vote. competing interests the authors declare that they have no competing interest. authors’ contributions ls carried out all analyses. pl helped with the r implementation and analysis. ls and sh developed the method and wrote the article. sh conceived of the study. all authors read and approved the final manuscript. supplementary material additional file 1 simulation study design. this file describes the simulation studies and presents r code used for simulating the data set. click here for file additional file 2 sensitivity and specificity of predictors in the disease gene expression data sets. for each data set and prediction method, the table reports the sensitivity and specificity estimated using 3-fold cross validation. more precisely, the table reports the average 3-fold cv estimate over random partitions of the data into folds. median sensitivity and specificity across data sets are summarized at the bottom. click here for file additional file 3 sensitivity and specificity of predictors in the uci machine learning benchmark data. for each data set and prediction method, the table reports the sensitivity and specificity estimated using 3-fold cross validation. more precisely, the table reports the average 3-fold cv estimate over random partitions of the data into folds. median sensitivity and specificity across data sets are summarized at the bottom. click here for file additional file 4 prediction accuracy when including pairwise interactions between features in the uci machine learning benchmark data. this table is an extension to table it shows the prediction accuracy of predictors other than rglm when considering pairwise interactions between features in the same uci mlbench data sets. although several predictors show improvement, none of them beats rglm.inter click here for file additional file 5 comparison of rglm based feature selection method with the rf based method of díaz-uriarte et al. for each data set in the disease gene expression data, the rf based variable selection method by díaz-uriarte et al selects a small set of genes. for each of the selected genes, the file reports the ranking in terms of the rglm variable importance measure timesselectedbyforwardregression. as expected, only a few of the selected genes have a high rank in terms of timesselectedbyforwardregression illustrating that these variable selection methods are different. click here for file additional file 6 effect of the number of bags on rglm predictor thinning. s this figure reports how prediction accuracy changes as variable thinning is applied to the rglm. results are averaged over the dichotomized gene traits in the mouse adipose data set. the five rows correspond to nbags values of respectively. within each row, the two panels have the same meaning as in figure click here for file additional file 7 prediction accuracy versus number of bags used for rglm. this figure presents the results for predicting gene traits in the brain cancer data set when different numbers of bags are used for constructing the rglm. each color represents one gene trait. binary outcome prediction. the gene traits were randomly selected from all gene traits used in the binary outcome prediction section. continuous outcome prediction. the gene traits were randomly selected from all gene traits used in the continuous outcome prediction. click here for file acknowledgements we acknowledge grant support from 1r01da030913- p50ca p30ca ul1tr we acknowledge the efforts of igc and expo in providing data set gse 
in contrast, a generalized linear model is very interpretable especially when forward feature selection is used to construct the model. this random generalized linear model predictor provides variable importance measures that can be used to define a “thinned” ensemble predictor that retains excellent predictive accuracy. rglm is a state of the art predictor that shares the advantages of a random forest with those of a forward selected generalized linear model. however, forward feature selection tends to overfit the data and leads to low predictive accuracy. a novel bootstrap aggregated glm predictor that incorporates several elements of randomness and instability often outperforms a host of alternative prediction methods including random forests and penalized regression models. comprehensive evaluations involving hundreds of genomic data sets, the uci machine learning benchmark data, and simulations are used to give glm based ensemble predictors a new and careful look. ensemble predictors such as the random forest are known to have superior accuracy but their black-box predictions are difficult to interpret. therefore, it remains an important research goal to combine the advantages of ensemble predictors with the advantages of forward regression modeling. to address this goal several articles have explored glm based ensemble predictors. since limited evaluations suggested that these ensemble predictors were less accurate than alternative predictors, they have found little attention in the literature.
 ensemble predictors such as the random forest are known to have superior accuracy but their black-box predictions are difficult to interpret. in contrast, a generalized linear model is very interpretable especially when forward feature selection is used to construct the model. however, forward feature selection tends to overfit the data and leads to low predictive accuracy. therefore, it remains an important research goal to combine the advantages of ensemble predictors with the advantages of forward regression modeling. to address this goal several articles have explored glm based ensemble predictors. since limited evaluations suggested that these ensemble predictors were less accurate than alternative predictors, they have found little attention in the literature.  comprehensive evaluations involving hundreds of genomic data sets, the uci machine learning benchmark data, and simulations are used to give glm based ensemble predictors a new and careful look. a novel bootstrap aggregated glm predictor that incorporates several elements of randomness and instability often outperforms a host of alternative prediction methods including random forests and penalized regression models. this random generalized linear model predictor provides variable importance measures that can be used to define a “thinned” ensemble predictor that retains excellent predictive accuracy.  rglm is a state of the art predictor that shares the advantages of a random forest with those of a forward selected generalized linear model. these methods are implemented in the freely available r software package randomglm. 