Abstract
Background
Genomewide association studies for complex diseases will produce genotypes on hundreds of thousands of single nucleotide polymorphisms (SNPs). A logical first approach to dealing with massive numbers of SNPs is to use some test to screen the SNPs, retaining only those that meet some criterion for futher study. For example, SNPs can be ranked by pvalue, and those with the lowest pvalues retained. When SNPs have large interaction effects but small marginal effects in a population, they are unlikely to be retained when univariate tests are used for screening. However, modelbased screens that prespecify interactions are impractical for data sets with thousands of SNPs. Random forest analysis is an alternative method that produces a single measure of importance for each predictor variable that takes into account interactions among variables without requiring model specification. Interactions increase the importance for the individual interacting variables, making them more likely to be given high importance relative to other variables. We test the performance of random forests as a screening procedure to identify small numbers of riskassociated SNPs from among large numbers of unassociated SNPs using complex disease models with up to 32 loci, incorporating both genetic heterogeneity and multilocus interaction.
Results
Keeping other factors constant, if risk SNPs interact, the random forest importance measure significantly outperforms the Fisher Exact test as a screening tool. As the number of interacting SNPs increases, the improvement in performance of random forest analysis relative to Fisher Exact test for screening also increases. Random forests perform similarly to the univariate Fisher Exact test as a screening tool when SNPs in the analysis do not interact.
Conclusions
In the context of largescale genetic association studies where unknown interactions exist among true riskassociated SNPs or SNPs and environmental covariates, screening SNPs using random forest analyses can significantly reduce the number of SNPs that need to be retained for further study compared to standard univariate screening methods.
Background
Genomewide association studies for complex diseases such as asthma, schizophrenia, diabetes, and hypertension will soon produce genotypes on hundreds of thousands of single nucleotide polymorphisms (SNPs). Due to the large number of SNPs tested and the potential for both genetic and environmental interactions, determining which SNPs modify the risk of disease is a methodological challenge. While the number of genotypes produced by candidate gene approaches will be somewhat less daunting, on the order of hundreds to thousands of SNPs, it will still be a considerable challenge to weed out the noise and identify the SNPs contributing to complex traits.
A logical first approach to dealing with massive numbers of SNPs is to first conduct univariate association tests on each individual SNP, in order to screenout those with no evidence for disease association. The primary goal of such a procedure is not to prove that a particular variant or set of variants influences disease risk, but to prioritize SNPs for further study. Using a univariate test at this stage will result in low power for SNPs with very small marginal effects in the population, even if the SNPs have large interaction effects. Of course, in addition to taking all individual SNPs, all SNP pairs could also be tested for association. However, when dealing with multiple thousands of SNPs at the outset, such an approach is cumbersome, and raises the question of where to stop: why not all sets of three, four, or even five SNPs as well?
Many modelbuilding methods exist for dealing with large numbers of predictors. For example, stochastic search variable selection (SSVS) [1], a form of Bayesian model selection, has been explored as a tool to discover joint effects of multiple loci in the context of genetic linkage studies [24]. However, these methods are limited in the number of predictors that can be included at one time, causing some researchers to resort to a twostage approach, in which only main effects are considered in a first stage, and interactions between loci with strong main effects are considered in a second stage. This approach can lead to the loss of important interactions with only weak main effects.
Multivariate adaptive regression splines (MARS) models have also been explored in the context of genetic linkage and association studies [5,6] with some degree of success. However, these and other model selection methods appear to be limited in the number of predictors that can reasonably be accommodated in one analysis, and the types of possible interactions that are allowed must be specified in advance. They are not suited to the initial task of identifying from a massive set of SNPs a subset for further analyses.
Combinatorial partitioning and multifactor dimensionality reduction [710] are closely related methods developed specifically to detect higherorder interactions among polymorphisms that predict trait variation. However, these methods are meant to identify interactions among small sets of SNPs, and have minimal power in the presence of genetic heterogeneity [10]. They are therefore inappropriate for use as a screening tool for searching through thousands of SNPs to identify those contributing to phenotypes in the context of wholegenome association studies. The problem remains: how do we reasonably weed down from thousands or hundreds of thousands of SNPs to a number that can be used by available modeling methods, without losing the interactions that we hope to model in the first place?
An additional concern to be considered is genetic heterogeneity. We define genetic heterogeneity to mean that there are multiple possible ways to acquire a disease or trait that can involve different subsets of genes. Traditional regression models are limited in their ability to deal with underlying genetic heterogeneity (see, e.g., [11]). If genetic heterogeneity also leads to phenotypic heterogeneity, then methods that classify individuals into phenotypic subgroups for further analysis can be successful. Likewise, if heterogeneity in genetic etiology is primarily due to ethnic background, subdividing samples by selfreported ethnicity or genetically defined subgroups can be a powerful antecedent to data analyses for the identification of complex disease genes. However, even in the realm of Mendelian genetic diseases, heterogeneity is rarely so simple. For example, multiple polymorphisms in each of two different genes are responsible for familial breast cancer in the relatively homogeneous subpopulation of Ashkenazi Jewish women [12]. When the root of the heterogeneity is not known a priori, traditional regression models, which lump all individuals into a single group and estimate average effects over the entire sample, are unlikely to successfully identify the genetic causes of diseases.
Classification trees and random forests
Treebased methods consist of nonparametric statistical approaches for conducting regression and classification analyses by recursive partitioning (see, e.g., Hastie et al. [13]). These methods can be very efficient at selecting from large numbers of predictor variables such as genetic polymorphisms those that best explain a phenotype. Tree methods are useful when predictors may be associated in some nonlinear fashion, as no implicit assumptions about the form of underlying relationships between the predictor variables and the response are made. They are welladapted to dealing with some types of genetic heterogeneity, as separate models are automatically fit to subsets of data defined by early splits in the tree.
The ease of interpretation of classification trees, along with their flexibility in accommodating large numbers of predictors and ability to handle heterogeneity, has resulted in increasing interest in their application to genetic association and linkage studies. Classification trees have been adapted for use with sibling pairs to subdivide pairs into more homogenous subgroups defined by nongenetic covariates [14], thus increasing the power to detect linkage under heterogeneity [15]. They have also shown promise for the dissection of complex traits for both linkage and association [16,17], and for exploring interactions [6]. A related adaptive regression method has also shown promise in selecting a small number of predictive SNPs from a set of hundreds of potential predictors [18]. Tree methods have also been used to identify homogeneous groups of cases for further analyses [19], and as an adjunct to more traditional association methods [20].
Classification trees are grown by recursively partitioning the observations into subgroups with a more homogeneous categorical response [21]. At each node, the explanatory variable (e.g., SNP) giving the most homogeneous subgroups is selected. Choosing alternative predictors that produce slightly suboptimal splits can result in very different trees that have similar prediction accuracy. The Random Forests methodology [22] builds on several other methods using multiple trees to increase prediction accuracy [2325]. A random forest is a collection of classification or regression trees with two features that distinguish it from trees built in a deterministic manner. First, the trees are grown on bootstrap samples of the observations. Second, a random selection of the potential predictors is used to determine the best split at each node. For each tree, a bootstrap sample is obtained by drawing a sample with replacement from the original sample of observations. The bootstrap sample has the same number of individuals as the original sample, but some individuals are represented multiple times, while others are left out. The leftout individuals, sometimes called "outofbag", are used to estimate prediction error. Because a different bootstrap sample is used to grow each tree, there is a different set of outofbag individuals for each tree. With a forest of classification trees, each tree predicts the class of an individual. For each individual, the predictions, or "votes", are counted across all trees for which the individual was outofbag, and the class with the most votes is the individual's predicted class. Random forests produce an importance score for each variable that measures its importance. This score can be used to prioritize the variables, much as pvalues from test statistics are used.
Using ensembles of trees built in this manner increases the probability that some trees will capture interactions among variables with no strong main effect. Unlike variable selection methods, interactions among predictors do not need to be explicitly specified in order to be utilized by a forest of trees. Instead, any interactions between variables serve to increase the importance of the individual interacting variables, making them more likely to be given high importance relative to other variables. Thus, random forests appear to be particularly wellsuited to address a primary problem posed by large scale association studies. In preliminary studies, we have shown the potential of random forests in the context of linkage analysis [26]. Other investigators are beginning to recognize the potential of the Random Forest methodology for studying SNP association [27] and classification [28].
To fully understand the basis of complex disease, it is important to identify the critical genetic factors involved, and to understand the complex relationships between genotypes, environment, and phenotypes. The few successes to date in identifying genes for complex disease suggest that despite carefully collected large samples, novel approaches are needed in the pursuit to dissect the multiple and varying factors that lead to complex human traits. Ultimately, the challenge in identifying polymorphisms that modulate the risk of complex disease is to find methods that can seamlessly handle large numbers of predictors, capitalize on and identify interactions, and tease apart the multiple heterogeneous etiologies. Here, we explore the use of the Random Forest methodology [22,29] as a screening tool for identifying SNPs associated with disease in the presence of interaction, heterogeneity, and large amounts of noise due to unassociated polymorphisms.
Results
Genetic models
We simulated complex diseases with sibling recurrence risk ratio for the disease (λ_{s}) fixed at 2.0 and population disease prevalence K_{p }equal to 0.10. These values are consistent with or lower than estimates from known complex genetic traits, such as Alzheimer disease, where estimates of cumulative prevalence in siblings of affected range from 30–40%, compared to a population prevalence of 10% at age 80 [30]. Such traits are understood to be caused by multiple interacting genetic and environmental factors. Our genetic models incorporate both genetic heterogeneity and multiplicative interaction as defined by Risch [31]: we simulate sets of 4, 8, 16, and 32 risk SNPs ("rSNPs") in linkage equilibrium, interacting in independent pairs or quartets to increase disease risk, and contributing equally to the overall sibling recurrence risk ratio of 2 and population disease prevalence of 0.10. For simplicity, we simulated the models such that each rSNP pair or quartet accounts for the same proportion of the genetic risk, and each SNP within a pair/quartet is responsible for an equal proportion of the genetic risk. Thus, all of the rSNPs simulated for a model have the same allele frequency and the same observed marginal effect in the population. We denote the models using the shorthand HhMm, where H = h (=2, 4, 8, 16) is the number of heterogeneous systems, and M = m (=2 or 4) is the number of multiplicatively interacting SNPs within each system. For example, 16 loci are responsible for the total λ_{s }= 2 and K_{p }= .10 for models H4M4 and H8M2, and 32 loci are responsible for models H8M4 and H16M2. Table 1 presents relevant features of our models. The Methods section describes the genetic models in more detail.
Table 1. Genetic models used for simulating casecontrol data.
Simulation and analysis
All analyses were performed on 100 replicate data sets of 500 cases and 500 controls. In addition to the rSNPs contributing to the trait, we simulated noise SNPs ("nSNPs"), independent of disease status, with allele frequencies distributed equally across the range .01–.99. To simulate the results of an association study, in which we do not expect to be lucky enough to genotype all polymorphisms related to a trait, we included only a subset of the total number of rSNPs in each analysis. We denote the analysis design using the shorthand KkSsNn, where K = k is the total number of rSNPs genotyped in the analysis, S = s is the number of SNPs within each interaction system genotyped, and N = n is the total number of SNPs genotyped in the design. Thus, NK is the total number of nSNPs in the analysis. For example, suppose the genetic model is H8M4, and the design is K4S2N100. Then out of the total of 8 × 4 = 32 rSNPs that contribute to the trait, four are genotyped: two interacting SNPs from within each of two heterogeneity systems. Six heterogeneity systems are not represented at all in the analysis. In addition to the four genotyped rSNPs, 1004 = 96 total nSNPs are also genotyped in the design.
Comparison of raw and standardized importance scores
Random forests version 5 software [29] produces both raw (I_{T}) and standardized (Z_{T}) variable importance scores (see Methods section for definitions of the scores). Little is known about the properties of importance indices under different distributions of the predictor variables. We use simulation to investigate their properties in the context of discrete predictors such as genetic polymorphisms conferring susceptibility to a complex trait.
We first compared the raw and standardized scores, in order to determine whether one might outperform the other in screening. We considered a K4S2N100 analysis design for each genetic model described in Table 1. I_{T }and Z_{T }are highly correlated; the average correlation coefficient over 100 replicate data sets ranged from .93 (H8M4) to >.99 (H2M2 and H4M2) (Table 2). The average correlation between the ranks based on I_{T }and Z_{T }for the 100 SNPs over the 100 replicate data sets was 0.98 for each of the six models (Table 2). Comparing the ranks of the four rSNPs among all SNPs, neither importance measure outperforms the other for all models (Table 3). The mean ranks of the rSNPs for the two measures are significantly different only for the H16M2 and H8M4 models. For H16M2, the average rank of the rSNPs is higher for Z_{T }than for I_{T}. The opposite is true for H8M4 (see Table 3).
Table 2. Summary of the correlation between I_{T }and Z_{T }("raw") and rank(I_{T}) and rank(Z_{T}) ("rank") for four rSNPs and 96 nSNPs over 100 replicate data sets: K4S2N100 analysis design.
Table 3. Comparison of ranks based on Z_{T }and I_{T }for the four rSNPs over 100 replicate data sets: K4S2N100 analysis design.
Ranking SNPs based on Z_{T }and Fisher pvalue
We next compared the ranking of rSNPs by importance score (Z_{T}) to ranking by Fisher Exact test pvalue using K4S2N100 and K4S2N1000 analysis designs, where two SNPs from each of the first two interaction systems are in the analysis. Figure 1 shows the proportion of replicates for which the top ranked 1, top 2, top 3, and top 4 SNPs are the four genotyped rSNPs in the data set for each of the four most complex genetic models. For N100, the random forest Z_{T }criterion ranks the four rSNPs as the most significant SNPs more often than the univariate Fisher Exact test association pvalue under all genetic models. The difference between the random forest and association pvalue ranking is less extreme for N1000. For the H8M4 genetic model, the results do not suggest that one ranking system is better than the other overall. Figure 2 shows the proportion of replicates for which all rSNPs are among the top N SNPs. In other words, it is the proportion of data sets for which none of the genotyped rSNPs are screened out, if the top ranking N SNPs are chosen for further study. For N100, a consistently higher proportion of replicates ranked using Z_{T }contain all of the rSNPs. Thus, for a given probability of retaining all of the rSNPs, more SNPs can be eliminated using the Z_{T }criterion than the Fisher exact test pvalue. For example, for model H16M2, only 15 SNPs must be retained to have 80% probability that the 4 rSNPs are in the retained set, while 44 SNPs must be retained if the pvalue criterion is used. The difference is less dramatic for H8M4: 37 SNPs give 80% probability that the four genotyped rSNPs are in the retained set for the Z_{T }criterion, compared to 43 for the pvalue criterion. For N1000, the advantage of the Z_{T }criterion is clear for the H8M2 and H16M2 models. For H4M4, the advantage of Z_{T }is minor, while for H8M4, ranking by Z_{T }appears to give poorer results than the pvalue criterion for the higher cutoff values of N. A second interpretation of Figure 2 is that, for any number of retained (not screenedout) SNPs, the probability that all of the genotyped rSNPs are retained is higher for the Z_{T }criterion than for the univariate pvalue criterion for all but the H8M4 model with 1000 total SNPs.
Figure 1. Proportion of replicates for which the most significant 1, 2, 3, and 4 SNPs are all rSNPs for K4S2N100 and K4S2N1000 analysis designs. Genetic models are listed on the plots. "RF" and "Fisher" refer to the random forest importance index Z_{T }and the Fisher Exact test pvalue. See text for notation description.
Figure 2. Proportion of replicates for which all rSNPs are among the topranking N SNPs for K4S2N100 and K4S2N1000 analysis designs. Other notation as in Figure 1.
Noticing that the analyses with all SNPs from an interacting system (e.g., the H16M2K4S2 simulations) had a more substantial improvement in ranking using Z_{T }over pvalue than the analyses with subsets of SNPs from an interacting system, we hypothesized that the interactions among the pairs of analyzed rSNPs influenced the improved ranking performance of the random forests over the univariate tests. To confirm this, we used the H8M4 genetic model and analyzed the data in the following manner. For a constant number of analyzed rSNPs included in the model (K = 4, 8, or 16) and a constant 96 nSNPs, we looked at the effect of increasing S, the number of rSNPs from each interaction system that were genotyped. Thus, for K8S1, along with 96 nSNPs, one SNP from each of the first 8 systems was included in the analysis, while for K8S4, all four SNPs in the first two systems were included in the analysis. For K8S3, three SNPs from the first two systems, and one from the third were included. Assuming that the random forest analysis was taking advantage of the interactions among the rSNPs, and that this was responsible for the improved performance of the random forests over the univariate tests, we expected the Fisher pvalues and random forest importance Z_{T }to perform similarly when only a single rSNP was genotyped from each system, and the random forests to perform increasingly better than the univariate Fisher tests as S increased from 1 to 4. Figures 3 and 4 show the results, which are consistent with this hypothesis. For the Fisher pvalues, the proportion of replicates for which the N most significant SNPs were rSNPs is similar for each S. For the random forest importance Z_{T}, the S = 1 analyses for each K were similar to the Fisher results, while for each increase in S, the proportion of replicates for which the N most significant SNPs were rSNPs increases (Figure 3) and the proportion of replicates for which all rSNPs are present at any cutoff point increases (Figure 4). The differences can be substantial: for the H8M4 model, with K = 4 rSNPs in the analysis, the number of most significant SNPs required to have 80% probability that all four rSNPs are included is 50, 34, 22, and 5, respectively for S1, S2, S2, and S4. We conclude that for a given number of rSNPs within a set of potential predictors, the more interacting SNPs there are, and the larger the groups of SNPs that interact, the greater the performance increase of the random forest analysis as compared to a univariate analysis.
Figure 3. Proportion of replicates for which the most significant N SNPs are all rSNPs. H8M4 genetic model. Analysis designs include 96 noise SNPs; K and S are listed on the plots. Other notation as in Figure 1.
Figure 4. Proportion of replicates for which all rSNPs are among the topranking N SNPs for H8M4 genetic model. Analysis designs include 96 noise SNPs; K and S are listed on the plots.
Magnitude of difference
Beyond simply ranking SNPs, we may wish to use the magnitude of the difference in importance or pvalue to determine which subset of topranked SNPs should be prioritized for further study. Thus, particularly for the cases where the rSNPs are among the topranked SNPs, we want to determine not just that Z_{T }ranks interacting rSNPs higher than the univariate test, but also that the differences in rank correspond to differences in magnitude of Z_{T }that are meaningful. In other words, we want to know how much "better" in terms of Z_{T }(or pvalue) the rSNPs are than the nSNPs. Toward this goal, we computed the difference between the importance Z_{T }of the top ranked rSNP and the top ranked nSNP:
D_{max}(Z_{T}) = max_{rSNP}(Z_{T})  max_{nSNP}(Z_{T}),
as well as the lowest ranked rSNP and the top ranked nSNP:
D_{min}(Z_{T}) = min_{rSNP}(Z_{T})  max_{nSNP}(Z_{T}).
Thus, D_{min}(Z_{T}) is positive when the lowest ranked rSNP is larger than the highest ranked nSNP, and negative when the lowest ranked rSNP is smaller than the highest ranked nSNP.
We computed the analogous quantities, D_{max}(log p) and D_{min}(log p), for the log_{10 }transformed Fisher Exact test pvalues. In Figure 5, we have plotted box plots of these differences for several models using analysis designs K4S2N100 and K4S2N1000. Pvalues for a paired Ttest of whether the mean difference is equal to 0 are also placed on the plot. For H8M2, the lowest ranking rSNP has Z_{T }that is significantly higher than the highest ranking nSNP for both N100 and N1000, while the difference in log10(p) is not significantly different from 0 for N100, and is significantly less than 0 for N1000. These plots illustrate that the positive differences are typically more extreme for Z_{T }than for log p, and that the negative differences are less extreme.
Figure 5. Distribution of difference in importance ZT between the top ranked rSNP and the top ranked nSNP (Dmax(ZT), and lowest ranked rSNP (Dmin(ZT)) and top ranked nSNP. Dmax(log p) and Dmin(log p): differences using log10 pvalue from the Fisher exact test. Beside each boxplot is the pvalue for the test of whether the mean difference over the 100 replicates is significantly different from 0. Genetic models listed in plot. Analysis design: K4S2, with N100 and N1000 shown on plot.
Discussion
A key advantage of the random forest approach is that the investigator does not have to propose a model of any kind. This is important in an initial genomewide or candidate region association study, where little is known about the genetic architecture of the trait. If interactions among SNPs exist, they will be exploited within the trees, and the variable importance scores will reflect the interactions. Therefore, we expect that when unknown interactions between true risk SNPs exist, the random forest approach to screening large numbers of SNPs will outperform a univariate ranking method in finding the risk SNPs among a large number of irrelevant SNPs. Our genetic models for simulation feature both multiplicative interaction and genetic heterogeneity. The multiplicative interaction results in a marginal effect in the population, the size of which is dependent on many factors, including the amount of heterogeneity. Thus, we have highlighted models in which univariate tests still have power, and shown that the random forest analysis can outperform these tests for selecting subsets of SNPs for further study. For models with genetic heterogeneity and interactions resulting in no main effect, similar to the models described by Ritchie et al. [10], the performance of random forests compares considerably more favorably to univariate tests (data not shown), since the univariate tests have no power when main effects are absent. Further investigation of how to determine a cutoff for SNPs to keep for further analysis is needed. Unfortunately, this task is likely to be strongly dependent on information that is impossible for an investigator to know a priori, such as the underlying genetic model and the ratio of associated risk SNPs to noise SNPs in an analysis.
Our results from analyses with four risk SNPs among 1000 SNPs suggest that even when a high proportion of the analyzed SNPs are unassociated, a random forest can rank interacting SNPs considerably higher than a univariate test, and that the proportional difference in importance between the risk SNPs and the best of the noise SNPs can be larger on average for a random forest. In our scaleup from 100 to 1000 total SNPs, we kept the number of risk SNPs constant. In practice, as we increase the number of SNPs genotyped, we expect that we will also increase the number of risk SNPs (or SNPs in linkage disequilibrium with risk SNPs) that are captured in an analysis. Thus, as a larger and large proportion of the genome or candidate region is captured by a scan, the more likely we will be to have all or most of sets of SNPs that interact, and thus the more likely we are to be in situations where random forest screening will outperform univariate screening of SNP data.
It is important to consider the tuning parameters for such analyses. Consistent with the recommendations made by Breiman and Cutler [22,29], the number of variables randomly selected at each split seems to have minimal effect over a wide range of values surrounding the square root of the number of covariates (SNPs). Breiman and Cutler do not recommend a method to determine the number of trees necessary for an analysis. The documentation examples typically use on the order of 100–1000 trees, but these examples are primarily in the context of prediction, without computing estimates of variable importance. In our experience with the simulated data sets presented here, in which the truly associated covariates are outnumbered considerably by those that are noise, multiple thousands of trees must be used in order to get stable estimates of the variable importance. In practice, we recommend building several forests for a data set with a given number of trees. If the ranking of variables by importance does not change significantly from forest to forest, then the number of trees is adequate.
We have examined the use of random forests in the context of association studies for complex disease with uncorrelated SNP predictors. Random forests can also be used when predictors are correlated, as is the case with multiple SNPs in linkage disequilibrium within a small genetic region. For any analysis procedure, the more highly correlated variables are, the more they can serve as surrogates for each other, weakening the evidence for association for any one of the correlated variables to the outcome. In a random forest analysis, limited simulations suggest that correlated variables lead to diminished variable importance for each correlated risk SNP (data not shown). One way to limit the problems presented by SNPs in linkage disequilibrium is to use haplotypes instead of SNPs as predictor variables in a random forest. Future challenges include quantifying more completely the effect of linkage disequilibrium among SNPs submitted to a random forest analysis, and developing random forests in the context of haplotypes.
Conclusions
With the increasing size of association studies, twostage analyses, in which in the first stage a subset of the loci are retained for further analyses, are becoming more common. The most frequently voiced concern for these analyses is that variables that interact to increase disease risk but have minimal main effects in the population will be missed. Random forest analyses address this concern by presenting a summary importance of each SNP that takes into account its interactions with other SNPs. Current implementations of random forests can accommodate up to one thousand of SNPs in one analysis with the computation of importance. Further, there is no reason to restrict the input variables to SNPs. Potential environmental covariates can also be easily accommodated, allowing for SNPs with no strong main effect, but environmental interactions, to be distinguished from unassociated SNPs. We have shown that when unknown interactions among SNPs exist in a data set consisting of hundreds to thousands of SNPs, random forest analysis can be significantly more efficient than standard univariate screening methods in ranking the true diseaseassociated SNPs highly. After identifying the topranked SNPs and other variables, and weeding out those unlikely to be associated with the phenotype, more thorough statistical analyses, including model building procedures, can be performed.
Methods
Variable importance
Rather than selecting variables for modeling, a random forest uses all available covariates to predict response. Here, we use measures of variable importance to determine which covariates (SNPs, in our case) or sets of covariates are important in the prediction. Breiman [22] proposed to quantify the importance of a predictor variable by disrupting the dependence between the variable and the response and measuring the change in the tree votes compared to the original observations. In practice, this is achieved by permuting the variable values among all observations in the outofbag sample of each tree. If the variable is predictive of the response, it will be present in a large proportion of trees and be near the root of those trees. Observations with a changed variable value may be directed to the wrong side of the tree, leading to vote changes from the right to the wrong class. Conversely, if the variable is not related to the response, it will be present in few trees and, when present, it will be near a terminal node, so that few tree votes will be changed. In Random Forests (version 5) Breiman and Cutler [29] define the importance index as follows. For individual i, let X_{i }represent the vector of predictor variable values, y_{i }its true class, V_{j}(X_{i}) the vote of tree j and t_{ij }an indicator taking value 1 when individual i is outofbag for tree j and 0 otherwise. Let X^{(A,j) }= (X_{1}^{(A,j)},..., X_{N}^{(A,j)}) represent the vector of predictor variables with the value of variable A randomly permuted among the outofbag observations for tree j, and X^{(A) }the collection of X^{(A,j) }for all trees where N is the total number of individuals in the sample. Letting 1(C) denote the indicator function taking value 1 when the condition C is true and 0 otherwise, the importance index averages over the trees of the forest, and is defined as:
where N_{j }represents the number of outofbag individuals for tree j and T is the total number of trees.
The importance index can be standardized by dividing by a standard error derived from the betweentree variance of the raw index I_{T}, . The standardized index is defined as:
The variance represents the tree to tree variance of I_{T}, rather than the variance of I_{T }due to the sampling of the individuals from a population: the magnitude of Z_{T }increases with the number of trees in the forest, and the number of trees is limited only by computing time. Thus, this standardized index cannot be treated as a Zscore in the traditional sense.
Simulation models and methods
For simplicity, assume each locus has the same effect, and let (q_{0}, q_{1}, q_{2}) represent the penetrance factors for 0, 1, or 2 risk alleles for an individual locus in a given model. Let G = {g_{11}, g_{12}, . ., g_{HM}} be the multilocus genotype for an individual, where g_{hm}(=0, 1, 2 risk alleles) indicates the individual's genotype at locus m (=1, . ., M) of heterogeneity system h (=1, . ., H). Then the penetrance for genotype G is defined as:
For example, for model H2M2, an individual with genotype G = {0101} would have penetrance P_{G }= 1  (1  q_{0}q_{1})^{2 }= 2q_{0}q_{1 } (q_{0}q_{1})^{2}.
The penetrance factors (q_{0}, q_{1}, q_{2}) and risk allele frequencies, as well as other features of our genetic models, are listed in Table 1. For a given model type, such as H4M4, and a given λ_{s }and K_{p}, there is a unique allele frequency when we make the assumption that each SNP subunit has equal effect (the given penetrance factor vector) in the population. We chose penetrance factors such that the risk alleles at each locus for the H2M2, H4M2, H8M2, and H16M2 models are approximately additive in effect on the penetrance factor scale. For H4M4 and H8M4, we chose penetrance factor vectors such that the risk alleles show a moderate degree of dominance.
The marginal genotype relative risks (GRRs) listed in Table 1 are the relative penetrances for heterozygote and homozygote carriers of each risk allele, as compared to noncarriers in the population, which would be observed if only a single rSNP were considered at a time. Thus, this is a measure of the observed effect size of each of the rSNPs in the population. The marginal GRRs are modest, in line with what might be expected when there are a large number of small effects contributing to a complex phenotype. For cases, the genotypes for pairs/quartets of SNPs within an interacting system are positively correlated, while SNPs from distinct systems are negatively correlated. The magnitude of the correlations decreases with increasing number of heterogeneity systems and increasing number of equaleffect SNPs interacting within each system.
Analysis
All analyses were performed on 100 replicate data sets of 500 cases and 500 controls. We treated the SNPs as ordinal predictors. Random forests have one primary tuning parameter: "mtry" the number of randomly picked covariates to choose among for each split. The Random Forest v5 manual [29] recommends trying the square root of the number of predictors, along with values smaller and larger than the square root, and choosing the value that minimizes the out of bag prediction error rate. We considered both the prediction error and the stability of the variable importance estimates when determining the values of mtry to use and the number of trees to grow. We found that the prediction error rate was very stable over a wide range of mtry for the number of trees we required for consistent measures of importance. We analyzed each replicate data set with 4–16 rSNPs and 96 nSNPs using a random forest of 5000 trees, choosing the best split from among a different randomlyselected set of 35 SNPs at each node (mtry = 35). On average, each replicate data set with 100 total SNPs took 40 minutes to complete on a 2.6 Ghz Intel Xeon processor. For data sets with 4 rSNPs and 996 nSNPs (1000 SNPs total), we used 15000 trees, and chose from among 125 SNPs at each node (mtry = 125). Analysis of each replicate of these data sets took 123 minutes on average. User time could potentially be substantially decreased by parallelprocessing: trees could be grown on separate nodes, and combined for analysis of importance. However, parallel treebuilding is not yet available in the Random Forest progream. To compare the performance of random forests with that of a univariate, oneSNPatatime approach, we tested for association between genotypes for each individual SNP and disease status using a Fisher Exact test [32].
Ranking of rSNPs
The random forest analysis produces the raw and standardized importance indices (I_{T}, Z_{T},), which can be used to rank order the importance of SNPs much as pvalues from association tests are used. Using either method, the ranking of the SNPs in an analysis can be used to prioritize which sets of SNPs (or genome regions) should be followed up with further genotyping and/or additional analyses. We use the convention that a rank of 1 is the highest ranking SNP.
We compare the ranking of the raw and standardized importance measures, and further compare these with the rank based on pvalues from a test of association, the Fisher Exact test, to determine whether the random forest can better discriminate susceptibility SNPs from SNPs unrelated to disease status when there is interaction and heterogeneity among and between SNPs.
Authors' contributions
KLL conceived of and led the design of the study, coordinated all phases of simulation and analysis, and drafted the manuscript. LBH and JS automated the simulation and random forest analysis procedures, and participated in the design and analysis of study results. PVE added critical insight to the development of the genetic models and participated in the interpretation of study results. All authors provided comments on a draft manuscript and read and approved the final manuscript.
Acknowledgments
We thank Kathleen Falls for her help with earlier versions of this work.
References

George EI, McCulloch RE: Variable Selection via Gibbs Sampling.
Journal of the American Statistical Association 1993, 88(423):881889.

Oh C, Ye KQ, He Q, Mendell NR: Locating disease genes using Bayesian variable selection with the HasemanElston method.
BMC Genet 2003, 4(Suppl 1):S69. PubMed Abstract  BioMed Central Full Text

Suh YJ, Ye KQ, Mendell NR: A method for evaluating the results of Bayesian model selection: application to linkage analyses of attributes determined by two or more genes.
Hum Hered 2003, 55:147152. PubMed Abstract  Publisher Full Text

Yi N, George V, Allison DB: Stochastic search variable selection for identifying multiple quantitative trait loci.
Genetics 2003, 164:11291138. PubMed Abstract  Publisher Full Text

York TP, Eaves LJ: Common disease analysis using Multivariate Adaptive Regression Splines (MARS): Genetic Analysis Workshop 12 simulated sequence data.
Genet Epidemiol 2001, 21 Suppl 1:S64954. PubMed Abstract

Cook NR, Zee RY, Ridker PM: Tree and spline based association analysis of genegene interaction models for ischemic stroke.
Stat Med 2004, 23(9):14391453. PubMed Abstract  Publisher Full Text

Nelson MR, Kardia SL, Ferrell RE, Sing CF: A combinatorial partitioning method to identify multilocus genotypic partitions that predict quantitative trait variation.
Genome Res 2001, 11(3):458470. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ritchie MD, Hahn LW, Roodi N, Bailey LR, Dupont WD, Parl FF, Moore JH: Multifactordimensionality reduction reveals highorder interactions among estrogenmetabolism genes in sporadic breast cancer.
Am J Hum Genet 2001, 69(1):138147. PubMed Abstract  Publisher Full Text

Hahn LW, Ritchie MD, Moore JH: Multifactor dimensionality reduction software for detecting genegene and geneenvironment interactions.
Bioinformatics 2003, 19(3):376382. PubMed Abstract  Publisher Full Text

Ritchie MD, Hahn LW, Moore JH: Power of multifactor dimensionality reduction for detecting genegene interactions in the presence of genotyping error, missing data, phenocopy, and genetic heterogeneity.
Genet Epidemiol 2003, 24(2):150157. PubMed Abstract  Publisher Full Text

Province MA, Shannon WD, Rao DC: Classification methods for confronting heterogeneity.
Adv Genet 2001, 42:273286. PubMed Abstract

LevyLahad E, Catane R, Eisenberg S, Kaufman B, Hornreich G, Lishinsky E, Shohat M, Weber BL, Beller U, Lahad A, Halle D: Founder BRCA1 and BRCA2 mutations in Ashkenazi Jews in Israel: frequency and differential penetrance in ovarian cancer and in breastovarian cancer families.
Am J Hum Genet 1997, 60(5):10591067. PubMed Abstract

Hastie T, Tibshirani R, Friedman JH: The elements of statistical learning : data mining, inference, and prediction. In Springer series in statistics. New York , Springer; 2001:xvi, 533.

Costello TJ, Swartz MD, Sabripour M, Gu X, Sharma R, Etzel CJ: Use of treebased models to identify subgroups and increase power to detect linkage to cardiovascular disease traits.
BMC Genet 2003, 4 Suppl 1:S66. PubMed Abstract  BioMed Central Full Text

Shannon WD, Province MA, Rao DC: Treebased recursive partitioning methods for subdividing sibpairs into relatively more homogeneous subgroups.
Genet Epidemiol 2001, 20(3):293306. PubMed Abstract  Publisher Full Text

Zhang H, Bonney G: Use of classification trees for association studies.
Genet Epidemiol 2000, 19(4):323332. PubMed Abstract  Publisher Full Text

Zhang H, Tsai CP, Yu CY, Bonney G: Treebased linkage and association analyses of asthma.
Genet Epidemiol 2001, 21 Suppl 1:S31722. PubMed Abstract

Kooperberg C, Ruczinski I, LeBlanc ML, Hsu L: Sequence analysis using logic regression.
Genet Epidemiol 2001, 21 Suppl 1:S62631. PubMed Abstract

Chang CJ, Fann CS: Using data mining to address heterogeneity in the Southampton data.
Genet Epidemiol 2001, 21 Suppl 1:S1805. PubMed Abstract

Wilcox MA, Smoller JW, Lunetta KL, Neuberg D: Using recursive partitioning for exploration and followup of linkage and association analyses.
Genet Epidemiol 1999, 17 Suppl 1:S3916. PubMed Abstract

Breiman L: Classification and regression trees. In The Wadsworth statistics/probability series. Belmont, Calif. , Wadsworth International Group; 1984:358 p..

Machine Learning 2001, 45:532. Publisher Full Text

Freund Y, Schapire R: Experiments with a new boosting algorithm.
Machine Learning: Proceedings of the Thirteenth International Conference 1996, 148156.

Schapire R, Freund Y, Bartlett P, Lee W: Boosting the margin: A new explanation for the effectiveness of voting methods.
The Annals of Statisics 1998, 26(5):16511686. Publisher Full Text

Bureau A, Dupuis J, Hayward B, Falls K, Van Eerdewegh P: Mapping complex traits using Random Forests.
BMC Genet 2003, 4 Suppl 1:S64. PubMed Abstract  BioMed Central Full Text

Horvath S, Kraft P: Using Random Forests to Detect Covariate Interaction Effects in CaseControl Studies: Applications to Screening for Disease Genes. In Joint Statistical Meetings. Volume Abstract #301988. San Francisco ; 2003.

Schwender H, Zucknick M, Ickstadt K, Bolt HM, The GENICA network: A pilot study on the application of statistical classification procedures to molecular epidemiological data.
Toxicol Lett 2004, 151.(1):291299. PubMed Abstract  Publisher Full Text

Breiman L, Cutler A: Random Forests. [http://statwww.berkeley.edu/users/breiman/RandomForests/cc_home.htm] webcite
Version 5 edition. 2004.

Farrer LA, Cupples LA: Determining the Genetic Component of a Disease. In Approaches to Gene Mapping in Complex Disease. Edited by L HJ, PericakVance MA. New York , John Wiley and Sons; 1998:93129.

Risch N: Linkage strategies for genetically complex traits. II. The power of affected relative pairs.
Am J Hum Genet 1990, 46(2):229241. PubMed Abstract

Agresti A: Categorical Data Analysis. 1st edition. New York , John Wiley & Sons; 1990:558.