The problem of locating multiple interacting quantitative trait loci (QTL) can be addressed as a multiple regression problem, with marker genotypes being the regressor variables. An important and difficult part in fitting such a regression model is the estimation of the QTL number and respective interactions. Among the many model selection criteria that can be used to estimate the number of regressor variables, none are used to estimate the number of interactions. Our simulations demonstrate that epistatic terms appearing in a model without the related main effects cause the standard model selection criteria to have a strong tendency to overestimate the number of interactions, and so the QTL number. With this as our motivation we investigate the behavior of the Schwarz Bayesian information criterion (BIC) by explaining the phenomenon of the overestimation and proposing a novel modification of BIC that allows the detection of main effects and pairwise interactions in a backcross population. Results of an extensive simulation study demonstrate that our modified version of BIC performs very well in practice. Our methodology can be extended to general populations and higher-order interactions.

---

POPULAR methods for mapping quantitative trait loci (QTL) include interval mapping (LANDER and BOTSTEIN 1989), composite interval mapping (ZENG 1993, 1994) and multiple QTL mapping (MQM; JANSEN 1993; JANSEN and STAM 1994). These statistical methods do not allow the location of QTL in situations when there are no main effects for the respective QTL, but there are (epistatic) interactions with other QTL (genes) that influence the quantitative trait. Epistatic QTL are known to play important roles in many disease studies, such as cancer (FIJNEMAN et al. 1996, 1998), and it is also suspected that they play a key role in the evolutionary process (WOLF et al. 2000).

A direct solution to detecting epistatic QTL is to search for several QTL simultaneously and fit an appropriate multiple regression model with interactions. However, the utility of such an approach, which is referred to as a multidimensional version of interval mapping, called multiple interval mapping (MIM; KAO et al. 1999), is limited by two interconnected issues. The first is the requirement of deciding how many terms (main effects and epistasis) should be included in the model. The second issue is the computational complexity of the search over the space of possible multidimensional models. To avoid these problems JANNINK and JANSEN (2001) and BOER et al. (2002) proposed one-dimensional genome searches as a means of mapping epistatic QTL. In particular they proposed an interesting extension of MQM by addressing a crucial problem pertaining to the choice of marker cofactors. By including all available markers in a regression equation and using a Bayesian approach to penalize large values of the corresponding regression coefficients many of the previously mentioned issues are eliminated. The disadvantage of this method is that, when detecting epistatic QTL, it requires the choice of "the effective dimension" (i.e., number of QTL) for epistatic interactions, which has strong influence on the power of detection.

An alternative way to approach the problem of mapping epistatic QTL relies on developing new methods for reducing the numerical complexity of MIM. In recent work CARLBORG et al. (2000), NAKAMICHI et al. (2001), and BROMAN and SPEED (2002) use random search methods to accelerate the search over the class of possible multidimensional models. The results from their approach hold great potential for further progress in solving the problem of the computational complexity of MIM.

Regardless of which method we use to search the genome for QTL we need to solve the problem of estimating QTL number, which in turn directly affects the dimensionality of the model space. The standard way of deciding how many main and interacting (QTL) effects should appear in the model relies on using many statistical tests (see KAO et al. 1999). A disadvantage of this approach is that it allows the comparison of only nested models. It is also unclear how to adjust the significance thresholds for each consecutive test.

Model selection criteria have been used as an alternative approach for the problem of model selection in QTL mapping. Two easy-to-compute model selection criteria that are often employed in statistics are the Akaike information criterion (AIC; AKAIKE 1974) or the Schwarz Bayesian information criterion (BIC; SCHWARZ 1978). These criteria belong to the family of the so-called penalized maximum likelihood methods and are based on the recommendation of choosing the model for which the likelihood of the data minus the penalty for the model dimension obtains the maximal value. These criteria were used by JANSEN (1993) and JANSEN and STAM (1994) to choose marker covariates for MQM and by PIEPHO and GAUCH (2001), NAKAMICHI et al. (2001), BALL (2001), BROMAN and SPEED (2002), and SIEGMUND (2003) to directly estimate QTL number. For a review and discussion of model selection methods as applied to QTL mapping see BALDING et al. (2002) or SILANPää and CORANDER (2002).

PIEPHO and GAUCH (2001) investigated many model selection criteria via simulation. In their study different criteria were used to choose pairs of markers flanking QTL. Their results suggest that out of the considered criteria BIC has the best properties and can be recommended for the estimation of the number of QTL with main effects. BROMAN and SPEED (2002), however, recommend a modification of BIC to select markers strongly associated with the trait. Contrary to PIEPHO and GAUCH (2001) they use BIC to choose single markers instead of pairs. BROMAN and SPEED (2002) observe that in this situation the original BIC has a tendency to overestimate the QTL number. To solve the problem of the overfitting they propose a modification of BIC, with a larger penalty for model dimension. Simulations reported in BROMAN and SPEED (2002) show that their modified version of BIC performs very well and detects the correct model more often than composite interval mapping does (ZENG 1993, 1994).

While both of the methods put forth by PIEPHO and GAUCH (2001) and BROMAN and SPEED (2002) can be used to estimate the number of QTL with main effects, they do not generalize directly to the situation where interaction terms appear in the model. Our extensive simulations (BOGDAN and DOERGE 2003) showed that the phenomenon of overfitting becomes even more significant when we allow interaction terms to appear in the model without the related main effects.

In the present work we concentrate on BIC, which, according to the QTL simulation study of PIEPHO and GAUCH (2001) and our independent simulations, performs better than other popularly used model selection criteria. In particular, we recall the Bayesian roots of BIC and explain the reasons why this criterion, when used to select single markers, has a tendency to overestimate the model dimension. To address this issue we follow the approach suggested by BALL (2001) and propose an easy modification of BIC that relies on taking into account the realistic prior distribution on the set of compared models. In comparison to BALL (2001) we extend the method to cover models with interactions and calibrate the prior to gain the control over the type I error of our procedure. An extensive simulation study verifies that our proposed criterion deals very well with the problem of overfitting the model and allows the detection of main effects and pairwise interactions in a backcross population. While our proposal is based on QTL mapping in a backcross population, our methodology can be extended to general populations and to higher-order interactions.

Consider a backcross population where q_ij denotes the genotype of the ith individual at the jth QTL: q_ij = –¹/₂ if the ith individual is homozygous at the jth QTL and q_ij = ¹/₂ if it is heterozygous. We assume that the relationship between the trait value Y_i and QTL genotypes is given by a normal regression model,

(1)

where m is the QTL number and _i N(0, ²) is the environmental noise. The second summation in our model corresponds to pairwise epistatic interactions. The formulation of the model allows some of the coefficients ß_j and _jl to be zero to accommodate cases when there are QTL that are not involved in epistatic effects. It also addresses the scenario when QTL might not have their own main effects, yet influence the quantitative trait by interacting with other genes, (i.e., epistasis). Later we use p to denote the number of QTL with main effects and q to denote the number of nonzero epistatic terms.

We rely on MIM (KAO et al. 1999) to simultaneously locate multiple QTL. This method requires fitting the model (1) for a dense grid of possible QTL positions. For each of the possible genomic locations the genotypes of the putative QTL are inferred using the genotypes of flanking markers and the EM algorithm (DEMPSTER et al. 1977) is employed to estimate parameters of the model (1). The locations for which the fitted model yields the largest likelihood are subsequently chosen.

A first step in the reduction of the complexity of MIM sometimes relies on identifying interesting genomic regions on the basis of an initial, relatively coarse search. In the Bayesian setting this approach was suggested by SEN and CHURCHILL (2001), who used an initial scan based on a 10-cM pseudo-marker grid. However, for the situation where an accurate genetic map exists a natural approach is to base the initial search on the net of marker positions and then use more refined methods (e.g., MIM) to search in the neighborhood of the chosen markers. BALL (2001), BROMAN and SPEED (2002), YI et al. (2003a), and XU (2003) successfully search over markers to locate multiple QTL and are justified in doing so on the basis of the fact that flanking markers absorb all the information associated with the QTL (WHITTAKER et al. 1996).

If we reduce MIM to a search over markers, then the problem of the QTL location reduces to the problem of choosing the best model of the form

(2)

where X_ij denotes the genotype of the ith individual at the jth marker; I is a certain subset of the set = {1, ... , N_m}, where N_m is the number of available markers; and U is a certain subset of x . For a backcross population the random variables X_iuX_iv correspond to the epistatic terms that are not correlated to any of the main effects. In particular, X_iuX_iv is not correlated to either X_iu or X_iv even if the uth and vth markers are statistically dependent via linkage. Thus, the epistatic effects are statistically not confounded with any of the main effects, and in most cases they will be detected only if the epistatic interactions are present.

One difficulty in fitting model (2) is the estimation of the number of main effects and interaction terms to be included in the model. There is a vast statistical literature on the choice of the number of terms in a linear model (see MILLER 1990 or MCQUARRIE and TSAI 1998) and there are many model selection criteria that can be used for this purpose. As mentioned earlier BROMAN and SPEED (2002) and PIEPHO and GAUCH (2001) recommend using the Schwarz BIC (SCHWARZ 1978) to estimate the number of QTL with main effects. In a general statistical context BIC recommends choosing the model that maximizes the expression

(3)

where is the vector of model parameters, L(Y|) is the likelihood of the data, k is the number of parameters (dimension of ), and n is the sample size. BIC belongs to the wide class of the so-called penalized maximum-likelihood methods and the second term in this criterion, ¹/₂k log n, is called the penalty for the complexity of the model. An important advantage of BIC is that for a wide range of statistical problems, and in particular for multiple regression, it is consistent (i.e., when the sample size grows to infinity, the probability of choosing the right model converges to 1). In the context of linear regression, maximizing S is equivalent to minimizing

(4)

where RSS is the residual sum of squares from regression.

Rationale for modifying BIC:

BROMAN and SPEED (2002) report that the original BIC, when used to select single markers with significant main effects, has a tendency to overestimate QTL number. On the basis of work not shown here (BOGDAN and DOERGE 2003) we have found that the tendency to overestimate QTL number becomes more significant when the portion (or entirety) of the genome under investigation increases. To understand this further we compare the rates at which the number of different models increases as the number of available markers increases. Our rationale is based on the observation that the number of possible models of the particular form (2), involving k distinct markers, is equal to , where N_m is the total number of available markers. Thus, when k is much smaller than N_m, the number of models involving k markers increases with N_m approximately like N_m^k. The difference in the numbers of possible "small" and "large" models increases quickly with N_m, and for large N_m the probability of choosing models with many components, just by random chance, is relatively high. Furthermore, for a large number of interaction terms, BOGDAN and DOERGE (2003) show that the original BIC has a tendency to choose models with epistatic terms even when in reality there is no epistasis.

The phenomenon of overestimation itself suggests the way the standard model selection criteria should be modified to make them useful for QTL mapping. Namely, the high rate at which the number of multidimensional models increases, when the number of available markers increases, suggests that the penalty for the model dimension should increase with this number. This condition is satisfied, for example, by criteria proposed by BROMAN and SPEED (2002) and SIEGMUND (2003). Second, the fact that there are many more interaction terms than the main effects suggests that the penalty for including an interaction should be larger than the penalty for including a main effect. Following these two suggestions we modify BIC by supplementing it with a realistic prior distribution on the set of possible models. Taking advantage of the fact that BIC is the approximation to the Bayesian rule for the choice of the "best" model we denote by _i = (µ, ß₁, ... , ß_p(i), ₁, ... , _q(i), ) the vector of parameters of the ith linear model, M_i, given by Equation 2. Here p(i) and q(i) denote the number of main effects and interaction terms involved in M_i. We assign a certain prior distribution for _i and denote the density of this distribution by f(_i). Moreover, let us denote the prior probability of the ith model by (i). Given that L(Y|_i, M_i) denotes the likelihood of the data given the vector of parameters _i, let p(Y|M_i) denote the likelihood of the data given the model M_i,

(5)

The posterior probability of the ith model, given the data, is

(6)

where l is the number of possible models.

The Bayesian rule recommends choosing the model for which the posterior probability P(M_i|Y) is the largest (see SCHWARZ 1978). Since the denominator in Equation 6 is the same for all considered models, Bayes' rule recommends choosing the model for which (i)p(Y|M_i) is the largest. The BIC criterion neglects the prior probabilities (i) of different models and approximates log p(Y|M_i) by log L(Y|_i, M_i) – ¹/₂(p(i) + q(i) + 2)log n, where _i is the maximum-likelihood estimator of _i, and p(i) + q(i) + 2 is the number of estimated parameters [i.e., p(i) + q(i) for main and epistatic effects, and 2 for µ and ]. Neglecting (i) corresponds to assigning the same prior probability to all considered models. While in many applications this approach is well justified, in the context of QTL mapping it lends itself to assigning unrealistically high prior probabilities to the events where many regressors are involved [e.g., when 200 markers are available, the number of different models involving 100 main effects is 9.05 x 10⁵⁸ and the prior probability of the event that 100 regressors are involved is >10⁵⁶ times larger than the prior probability of the event that there is just one regressor]. Motivated to improve on this we suggest supplementing BIC with a more realistic prior distribution, , on the class of possible models, and choosing the model for which

(7)

obtains a maximum.

In the context of multiple regression

where C(n) is the constant dependent only on n, and maximizing (7) is equivalent to minimizing the quantity

Prior distribution :

Assume N_m markers are available, and therefore N_m potential regressors and N_e = (N_m(N_m – 1))/2 potential interaction terms. The number of all models of the form (2) that can be constructed using subsets of N_m markers is equal to 2^N_m+N_e. To assign prior probabilities to these models we follow the standard solution proposed in GEORGE and MCCULLOCH (1993). Namely, we assign the probability to the event that the ith main effect appears in the model and probability to the event that the jth interaction term appears in the model. Our prior distribution assumes that particular terms enter the model independently of others and for a particular model M_i involving p(i) main effects and q(i) interactions we obtain

This choice of prior implies that the prior distributions on the number of main effects and epistatic terms are binomial with parameters N_m and , and N_e and , respectively.

For simplicity we consider and as = 1/l, = 1/u, where l and u are certain natural numbers, and restate the prior distribution as

where C(N_m, N_e, l, u) is a constant dependent on N_m, N_e, l, and u. Incorporating this prior distribution into the BIC [modified Schwarz BIC (mBIC)] allows the following rule: choose the model that minimizes

(8)

The expected values of the prior distribution for the number of main effects are equal to N_m/l and N_e/u for the number of interaction terms. Therefore, since the choice of l and u should reflect our prior knowledge on the QTL number, the values of l and u should be relatively small when we expect many QTL and large when we expect only few. Extensive simulations were performed for the purpose of investigating the standard values of l and u when we have no prior knowledge on the QTL number. We let l and u take on values in such a way that for the sample sizes n 200 the probability of type I error (detecting at least one QTL when there are none) does not exceed 5%. We observed that when markers are densely spaced (distance between markers is not >20 cM) we can obtain our aim by keeping the expected values of the number of main effects and interaction terms at a constant level close to 2. In particular, and as is seen next, in our simulations we used values l N_m/2.2 and u N_e/2.2. In the APPENDIX we present results of some theoretical calculations that support our empirical choice of l and u. These calculations yield approximate bounds on the type I error of our procedure and demonstrate that the proposed choice of l and u solves the problem of multiple comparisons and allows control of the type I error. In comparison to the original BIC the penalty in our proposed/modified criterion involves additional terms 2p(i)log((N_m/2.2) – 1) and 2q(i)log((N_e/2.2) – 1). A similar additional penalty appears in the criterion proposed by SIEGMUND (2003), who approaches the problem of QTL mapping differently by treating it as a change-point problem. These additional terms make our criterion similar to the risk inflation criterion (RIC) proposed by FOSTER and GEORGE (1994) in which the penalty for including k orthogonal regressors is equal to 2k log t, where t is the total number of available regressors. Note, however, that when n tends to infinity these additional terms are overshadowed by the BIC penalty (p(i) + q(i)) log n and, contrary to RIC, our criterion has the asymptotic properties of the BIC (i.e., consistency).

We employ computer simulations to evaluate the applicability of our proposed modification to the BIC criterion. Marker and QTL genotypes are simulated for a backcross population using 12 chromosomes of the length 100 cM for sample sizes n = 200 and n = 500. The number of QTL with main effects ranges between 0 and 12, and the number of epistatic terms between 0 and 5 (Tables 1 and 2) . Models 4, 5, and 11–14 (Table 1) are included to allow for a direct comparison to the results of BROMAN and SPEED (2002), as indicated by model 12, and to the results of PIEPHO and GAUCH (2001)(models 4, 5, 11, 13, and 14). Since we are interested in how our proposed criterion adjusts to the number of available markers, we search for QTL over 1, 5, and 12 chromosomes and use marker spacings of 5, 10, and 20 cM. The number of available markers and interaction terms, as well as the corresponding values of l and u for each of these experiments, is specified in Table 3 . The forward selection procedure (see, e.g., MILLER 1990) is used to search the space of possible multidimensional models. At each consecutive step we test all terms (main and interaction) not yet in the model and choose the one whose presence in the model yields the lowest value of the modified BIC criterion (Equation 8; mBIC). To save computational time the procedure is stopped after 30 steps and the resulting 31 models are evaluated on the basis of minimizing the mBIC (Equation 8). The number of steps is restricted to 30 since the largest model we use in the simulations has only 12 terms. Actually, we observe that for all the cases that we considered, the mBIC criterion was minimized by models with <20 terms and that increasing the number of steps above 20 had no influence on the results. However, in real data studies, when one does not want to bound the QTL number, we suggest using a larger number of steps.

View this table: In this window In a new window   TABLE 1 Simulation models

 

View this table: In this window In a new window   TABLE 2 Details of epistatic effects employed in simulation (Table 1)

 

View this table: In this window In a new window   TABLE 3 Penalty coefficients l and u used in the modified BIC (mBIC)

 

The results of searching over 1, 5, and 12 100-cM chromosomes, respectively, with markers spaced every 10 cM are shown in Tables 4– 6 , while Table 7 reports the results for varying marker distances. The number of correctly identified terms (corr. id.), averaged across 100 simulations, and the average number of false positives (extr.) are reported. The false positives that occur are divided into categories depending on their linkage to true QTL. Following PIEPHO and GAUCH (2001) we classify the main effect to be correct if it corresponds to a marker lying within 15 cM of the true QTL. If two markers from the neighborhood of one QTL are chosen, one of these markers is arbitrarily classified as extraneous. Epistatic terms are classified as correct if both markers involved lie within 15 cM of the true QTL. For the no-QTL model (1) the percentage of replicates for which the model with no QTL was chosen is reported. While the 15-cM margin is somewhat arbitrary it accommodates our situation well and illustrates the performance of our criterion. Recall that our main goal is the estimation of QTL number and not the precise location of QTL. If we use a narrower range (i.e., <15 cM), then some of the properly identified terms will be classified as extraneous due to the relatively large error of localization of weak QTL that is inherent to all QTL mapping procedures.

View this table: In this window In a new window   TABLE 4 Results from 100 simulations that each search over 12 100-cM chromosomes with markers spaced every 10 cM

 

View this table: In this window In a new window   TABLE 6 Results from 100 simulations that each search over five 100-cM chromosomes with markers spaced every 10 cM

 

View this table: In this window In a new window   TABLE 5 Results from 100 simulations that each search over one 100-cM chromosome with markers spaced every 10 cM

 

View this table: In this window In a new window   TABLE 7 Results of the search over 12 100-cM chromosomes based on 100 simulations and the sample size n = 200

 
Our modification to BIC performs very well (Tables 4–7) in practice, adjusts appropriately to the number of available markers under consideration, and rarely overestimates. Furthermore, in all of the examples we considered the probability of incorrectly detecting at least one QTL, when there are none, does not exceed 0.06 and the average number of extraneous QTL, which are not linked to true QTL, rarely exceeds 0.10. We also observe that the average number of extraneous epistatic terms never exceeded 0.05. This confirms our expectations that in the backcross population epistatic effects are usually detected only when they really exist. Since we set the expected values of the prior distribution for the number of main effects and interaction terms to be equal to 2.2, our criterion more easily identifies models with a small number of terms. The properties of our proposed criterion quickly improve with increasing sample size. Therefore, the accuracy of detecting small models increases (see models 1, 6, and 7 in Table 4) as does the ability to correctly identify models with larger numbers of QTL (see models 12, 13, 16, and 17 in Table 4). We are aware that the chance of correctly identifying QTL depends on its heritability. In other words, when the variance of the error is equal to 1.0 and the sample size is n = 200, our criterion usually detects main effects with coefficients ß 0.76 (the heritability of the single QTL with such a ß is 0.13) and interaction terms with 2 (broad sense heritability of 0.20 with just one such epistatic term in the model) even when they appear in larger models. When the sample size is increased to n = 500 our criterion usually detects main effects with ß 0.50 (individual h² 0.06) and interaction terms with 1.5 (individual h² 0.12). The proposed criterion (mBIC) works particularly well if QTL are located close to markers (compare models 4 and 6, and 5 and 7, in Tables 4– 6 and models 4 and 8 in Table 7). When QTL are located in the middle of an interval defined by two markers it is sometimes the case that both flanking markers are chosen, which partially explains the relatively large number of false positives for models 4 and 15. An additional reason for the sometimes larger number of extraneous linked QTL is a statistical error of localization of weak QTL. In some cases the correct model was appropriately identified, but the chosen markers were slightly farther apart from the true QTL than our set limit of 15 cM. On the basis of this reasoning some of the false positives correspond to correctly identified, but incorrectly localized QTL. Comparing results of our simulations with the results reported in PIEPHO and GAUCH (2001) and BROMAN and SPEED (2002) we observe, for models with only main effects, that our modification of BIC (mBIC) performs similarly to the criteria proposed in these earlier articles. More importantly, however, our criterion allows the detection of epistatic terms whereas the criteria of PIEPHO and GAUCH (2001) and BROMAN and SPEED (2002) do not.

The method proposed in this article can be viewed as a simplification of standard Bayesian methods used for QTL mapping. In a series of articles SATAGOPAN and YANDELL (1996), SATAGOPAN et al. (1996), HEATH (1997), UIMARI and HOESCHELE (1997), SILANPää and ARJAS (1998), STEPHENS and FISCH (1998), and YI and XU (2000) use the full Bayesian approach and Markov chain Monte Carlo simulations to estimate posterior distributions of QTL locations and other parameters in the regression model. YI et al. (2003a), XU (2003), and KILPIKARI and SILANPää (2003) reduce the number of parameters generated by Markov chain Monte Carlo (MCMC) by restricting the search to marker positions. YI and XU (2002) and YI et al. (2003b) extend the standard Bayesian MCMC approach to search for epistatic QTL. The common feature shared by the works of these authors is that they require multiple generations from the conditional distributions of all parameters in the regression model and are very computationally demanding. Moreover, as noted by BALL (2001), "a major challenge remains to obtain a rapidly converging sampler for the full Bayesian model." SEN and CHURCHILL (2001) avoided using MCMC by employing an independent sample Monte Carlo approach to generate multiple versions of pseudo-marker genotypes on the dense grid of genomic locations. They computed weights for each pseudo-marker realization by integrating out parameters of the related regression models and then used them to approximate the posterior distribution of the QTL locations. Our method, similar to the methods of BALL (2001) and BROMAN and SPEED (2002), is a further simplification of Bayesian methodology and seems to be particularly useful when one needs to search over a large space of possible models with interactions.

The modified BIC that is presented here is closer than the original BIC to the concept of Bayesian thinking since it introduces the prior distribution on the number of main effects and epistatic terms. We concentrate mainly on the situation when there are no specific expectations on the number of QTL and calibrate the prior so as to gain control over the type I error of our procedure. However, we strongly suggest that in the case when some prior information is available it should be included and the penalty should be adjusted accordingly. To estimate the type I error in that case one could use computer simulations or the permutation method of CHURCHILL and DOERGE (1994).

In principle, the modified version of BIC suggested in this article could be used to approximate posterior probabilities of different models according to the formula

(9)

where l is the number of possible models (see also BALL 2001). While we are very much aware of the importance of this formulation, which could allow one to estimate the uncertainty related to the choice of the best model and to use Bayesian averaging to estimate main and epistatic effects, we point out that due to the huge number of possible models with interactions it is practically impossible to compute its denominator. To reduce the number of terms in the Equation 9 one could apply Occam's window algorithm proposed by RAFTERY et al. (1997), which relies on discarding models that receive little support from the data. However, the corresponding search procedure proposed in MADIGAN and RAFTERY (1994) seems to be inadequate in our setting due to the large number of nonnested models. In practice one may reduce the number of models considered by performing a separate search for each pair of chromosomes, which in turn is usually good enough to detect pairwise interactions. But even in this case, the number of possible models with interactions will usually be too large to apply Equation 9.

To solve the problem of multiplicity of models and to identify the best one, we applied forward selection procedure, which is simple and quick. Our simulations, as well as results reported in BROMAN and SPEED (2002), show that forward selection performs very well in this setting. We are, however, aware that there are some particular cases (and a real analysis is always a particular case) when the forward selection procedure does not detect the best model. Thus, although statistically we do not expect much improvement by replacing forward selection with a more refined search strategy, we still recognize the need for further research in this direction.

Although this article is concerned solely with detecting main effects and pairwise interactions, theoretically the proposed method can be directly generalized to identify higher-order interactions. To retain control over the type I error of the corresponding procedure, it is anticipated that higher-order interactions should be penalized even more than pairwise epistatic terms. However, the utility of this approach needs to be verified by additional research, since there are two main difficulties related to any extensions of our work. First is the numerical complexity of the search over a rapidly increasing number of models with higher-order interactions, which can most likely be addressed by developing a suitable search strategy and increasing computer power. The second issue is more difficult and of a more theoretical nature. If we do not have prior expectations on the number of main and epistatic effects the method outlined in this article can be used to control the overall type I error. In this case, when we increase the potential number of regressors by including higher-order interactions, we must also increase the penalties for main effects and pairwise interactions. Thus, an attempt to detect higher-order interactions will result in decreasing power of detection of simpler effects and can be offset only by larger sample sizes. When some prior information on the number of main effects and interactions is available the power will be less affected since the method can be used in a subjective way via an appropriate adjustment of the penalties.

In this article we did not address the problem of missing marker data. Currently in the QTL mapping literature three methods exist, which are designed to solve this problem by using genotypes of neighboring markers. They include HALEY and KNOTT (1992) regression, the E-M algorithm of JANSEN and STAM (1994), or multiple imputations of missing genotypes proposed by SEN and CHURCHILL (2001) and BALL (2001). We believe that the application of any of these methods will leave the mBIC unaffected by a moderate proportion of missing marker data. The missing data methods can also be used to apply mBIC to search for QTL within intermarker intervals.

The method proposed in this article selects markers strongly associated with the trait and does not explicitly use the information from the distance between them. Therefore, in principle the mBIC approach is not sensitive to map errors. However, the application of any of the missing data methods will make our method sensitive to map errors in the same way as standard interval mapping. Our method can be also influenced by selective genotyping and genotyping errors, since selective genotyping will change the correlation structure in the design matrix and might result in partial confounding of epistatic and main effects. However, our approach is able to select the proper markers out of many strongly correlated neighbors; therefore we believe that it is also robust to any partial confounding of main and epistatic effects. The influence of genotyping errors will depend on the marker information that is affected. In our mBIC criterion, as well as in other standard model selection criteria, the information on the data appears only in RSS. Thus, we do not expect a significant difference between our criterion and others with respect to the sensitivity to genotyping errors.

Locating QTL, and more importantly their interactions, remains an open problem in both the QTL mapping and statistical communities. Current multiple interval mapping methods are plagued by two intimately related issues. First is the problem of estimating the number of QTL and their interactions. And second is the related issue of searching over the space of all possible multidimensional models that comprise the computationally complex space. Realizing that the second issue is impacted by the approach of the first issue, we have presented a model selection criterion that allows more accurate assessment than the original BIC criterion, from which we started. Using simulations in a backcross setting we demonstrate that the mBIC does well to locate multiple interacting QTL. Extensions of our method to more general designs are possible and are currently under investigation. Furthermore, the proposed criterion can be used outside the context of genetics to estimate the number of additive effects and interactions in the general framework of multiple regression.

BOUND FOR THE TYPE I ERROR

Our procedure recommends choosing the model that maximizes the criterion

(A1)

The number of all possible one-dimensional models [models for which p(i) + q(i) = 1] is equal to N_m + N_e, where, as before, N_m is the number of available markers and N_e = (N_m(N_m – 1))/2 is the number of possible interactions. Let ₁ denote the maximum of the criterion (A1) over all such one-dimensional models and let ₀ = log L₀(Y|, ) be the value of the criterion for the null model involving no markers (p + q = 0). Let D = p + q be the number of terms in the model chosen by our procedure. It holds that

We bound the probability of the first, dominating term of the right-hand side of the above equality, under the null hypothesis of no QTL.

Consider a given one-dimensional model M_i and a corresponding value of our criterion

The model M_i will be preferred over the model with no QTL if _Mi > ₀, or equivalently

(A2)

Under the null hypothesis of no QTL 2 logL(Y|_i)/L₀(Y|, ) has asymptotically ² distribution with 1 d.f. (SERFLING 1980). Thus, P is asymptotically equivalent to

where Z is a N(0, 1) random variable.

Note that ₁ > ₀ if Equation A2 holds for at least one of the one-dimensional models. Therefore, by Bonferroni inequality, for any > 0 and sufficiently large n it holds

(A3)

For each x > 0 it holds that

Thus (A3) yields

For the proposed values l = N_m/2.2 and u = N_e/2.2 the right-hand side of the above inequality is approximately equal to

(A4)

Using the proposed values of l and u allows one to eliminate N_m and N_e from the bound numerator and thus helps solve the multiple-comparisons problem.

For the sample size n = 200 and N_m and N_e used in our experiments the bound given by (A4) takes values from the interval between 0.0574 (for N_m = 252 and N_e = 31,626; 12 100-cM chromosomes with markers spaced every 5 cM) and 0.0801 (for N_m = 11 and N_e = 55; one 100-cM chromosome with markers spaced every 10 cM), which gives a satisfactory approximation for the empirical type I error obtained from simulations.

We thank Jim Berger, Andreas Futschik, and Joanna Szyda for helpful discussions and two anonymous reviewers for comments that greatly improved the presentation of the results. We also thank Adam Zagdanski and Monika Horobiowska-Kaczmarz for help in performing simulations. The results reported in this article were presented at the Seventh Purdue Symposium on Statistics, June 19–24, 2003. R.W.D. is funded in part by grants from the U.S. Department of Agriculture-Initiative for Future Food and Agricultural Systems (00-52100-9615) and the National Science Foundation (0115109-MCB).

AKAIKE, H., 1974 A new look at the statistical model identification. IEEE Trans. Automat. Control AC-19: 716–723.[CrossRef]

BALDING, D., A. D. CAROTHERS, Y. L. MARCHINI, L. R. CARDON, A. VETTA et al., 2002 Discussion on the meeting on statistical modelling and analysis of genetic data. J. R. Stat. Soc. B 64: 737–775.[CrossRef]

BALL, R., 2001 Bayesian methods for quantitative trait loci mapping based on model selection: approximate analysis using the Bayesian information criterion. Genetics 159: 1351–1364.[Abstract/Free Full Text]

BOER, M. P., C. J. F. TER BRAAK and R. C. JANSEN, 2002 A penalized likelihood method for mapping epistatic quantitative trait loci with one-dimensional genome searches. Genetics 162: 951–960.[Abstract/Free Full Text]

BOGDAN, M., and R. W. DOERGE, 2003 Mapping multiple interacting quantitative trait loci with multidimensional genome searches. Technical Report 04-03. Department of Statistics, Purdue University, West Lafayette, IN.

BROMAN, K. W., and T. P. SPEED, 2002 A model selection approach for the identification of quantitative trait loci in experimental crosses. J. R. Stat. Soc. B 64: 641–656.[CrossRef]

CARLBORG, Ö., L. ANDERSSON and B. KINGHORN, 2000 The use of a genetic algorithm for simultaneous mapping of multiple interacting quantitative trait loci. Genetics 155: 2003–2010.[Abstract/Free Full Text]

CHURCHILL, G. A., and R. W. DOERGE, 1994 Empirical threshold values for quantitative trait mapping. Genetics 138: 963–971.[Abstract]

DEMPSTER, A. P., N. M. LAIRD and D. B. RUBIN, 1977 Maximum likelihood from incomplete data via EM algorithm. J. R. Stat. Soc. B 39: 1–38.

FIJNEMAN, R. J. A., S. S. DE VRIES, R. C. JANSEN and P. DEMANT, 1996 Complex interactions of new quantitative trait loci, Sluc1, Sluc2, Sluc3, and Sluc4, that influence the susceptibility to lung cancer in the mouse. Nat. Genet. 14: 465–467.[CrossRef][Medline]

FIJNEMAN, R. J. A., R. C. JANSEN, M. A. VAN DER VALK and P. DEMANT, 1998 High frequency of interactions between lung cancer susceptibility genes in the mouse: mapping of Sluc5 to Sluc14. Cancer Res. 58: 4794–4798.[Abstract/Free Full Text]

FOSTER, D. P., and E. I. GEORGE, 1994 The risk inflation criterion for multiple regression. Ann. Stat. 22: 1947–1975.

GEORGE, E. I., and R. E. MCCULLOCH, 1993 Variable selection via Gibbs sampling. J. Am. Stat. Assoc. 88: 881–889.[CrossRef]

HALEY, C. S., and S. A. KNOTT, 1992 A simple regression method for mapping quantitative trait loci in line crosses using flanking markers. Heredity 69: 315–324.[Medline]

HEATH, S. C., 1997 Markov chain Monte Carlo segregation and linkage analysis for oligogenic models. Am. J. Hum. Genet. 61: 748–760.[Medline]

JANNINK, J.-L., and R. JANSEN, 2001 Mapping epistatic quantitative trait loci with one-dimensional genome searches. Genetics 157: 445–454.[Abstract/Free Full Text]

JANSEN, R. C., 1993 Interval mapping of multiple quantitative trait loci. Genetics 135: 205–211.[Abstract]

JANSEN, R. C., and P. STAM, 1994 High resolution of quantitative traits into multiple loci via interval mapping. Genetics 136: 1447–1455.[Abstract]

KAO, C.-H., Z-B. ZENG and R. D. TEASDALE, 1999 Multiple interval mapping for quantitative trait loci. Genetics 152: 1203–1216.[Abstract/Free Full Text]

KILPIKARI, R., and M. J. SILANPää, 2003 Bayesian analysis of multilocus association in quantitative and qualitative traits. Genet. Epidemiol. 25: 122–135.[CrossRef][Medline]

LANDER, E. S., and D. BOTSTEIN, 1989 Mapping Mendelian factors underlying quantitative traits using RFLP linkage maps. Genetics 121: 185–199.[Abstract/Free Full Text]

MADIGAN, D., and A. E. RAFTERY, 1994 Model selection and accounting for model uncertainty in graphical models using Occam's window. J. Am. Stat. Assoc. 89: 1535–1546.[CrossRef]

MCQUARRIE, A. D. R., and C.-L. TSAI, 1998 Regression and Time Series Model Selection. World Scientific Publishers, Singapore.

MILLER, A. J., 1990 Subset Selection in Regression. Chapman & Hall, London.

NAKAMICHI, R., Y. UKAI and H. KISHINO, 2001 Detection of closely linked multiple quantitative trait loci using genetic algorithm. Genetics 158: 463–475.[Abstract/Free Full Text]

PIEPHO, H.-P., and H. G. GAUCH, JR., 2001 Marker pair selection for mapping quantitative trait loci. Genetics 157: 433–444.[Abstract/Free Full Text]

RAFTERY, A. E., D. MADIGAN and J. A. HOETING, 1997 Bayesian model averaging for linear regression models. J. Am. Stat. Assoc. 92: 179–191.[CrossRef]

SATAGOPAN, J. M., and B. S. YANDELL, 1996 Estimating the number of quantitative trait loci via Bayesian model determination. Special Contributed Paper Session on Genetic Analysis of Quantitative Traits and Complex Diseases, Biometric Section, Joint Statistical Meetings, Chicago.

SATAGOPAN, J. M., B. S. YANDELL, M. A. NEWTON and T. C. OSBORN, 1996 Bayesian model determination for quantitative trait loci. Genetics 144: 805–816.[Abstract]

SCHWARZ, G., 1978 Estimating the dimension of a model. Ann. Stat. 6: 461–464.

SEN, S., and G. A. CHURCHILL, 2001 A statistical framework for quantitative trait mapping. Genetics 159: 371–387.[Abstract/Free Full Text]

SERFLING, R., 1980 Approximation Theorems of Mathematical Statistics. Wiley, New York.

SIEGMUND, D., 2003 Model selection in irregular problems: applications to mapping QTLs. Biometrika (in press).

STEPHENS, D. A., and R. D. FISCH, 1998 Bayesian analysis of quantitative trait locus data using reversible jump Markov chain Monte Carlo. Biometrics 54: 1334–1347.[CrossRef]

SILANPää, M. J., and E. ARJAS, 1998 Bayesian mapping of multiple quantitative trait loci from incomplete inbred line cross data. Genetics 148: 1373–1388.[Abstract/Free Full Text]

SILANPää, M. J., and J. CORANDER, 2002 Model choice in gene mapping: what and why. Trends Genet. 18: 301–307.[CrossRef][Medline]

UIMARI, P., and I. HOESCHELE, 1997 Mapping linked quantitative trait loci using Bayesian analysis and Markov chain Monte Carlo algorithms. Genetics 146: 735–743.[Abstract]

WHITTAKER, J. C., R. THOMPSON and P. M. VISSCHER, 1996 On the mapping of QTL by regression of phenotype on marker-type. Heredity 77: 23–32.

WOLF, J. B., E. D. BRODIE, III and M. J. WADE (Editors), 2000 Epistasis and the Evolutionary Process. Oxford University Press, New York.

XU, S., 2003 Estimating polygenic effects using markers of the entire genome. Genetics 163: 789–801.[Abstract/Free Full Text]

YI, N., and S. XU, 2000 Bayesian mapping of quantitative trait loci for complex binary traits. Genetics 155: 1391–1403.[Abstract/Free Full Text]

YI, N., and S. XU, 2002 Mapping quantitative trait loci with epistatic effects. Genet. Res. 79: 185–198.[CrossRef][Medline]

YI, N., V. GEORGE and D. B. ALLISON, 2003a Stochastic search variable selection for identifying multiple quantitative trait loci. Genetics 164: 1129–1138.[Abstract/Free Full Text]

YI, N., S. XU and D. B. ALLISON, 2003b Bayesian model choice and search strategies for mapping interacting quantitative trait loci. Genetics 165: 867–883.[Abstract/Free Full Text]

ZENG, Z-B., 1993 Theoretical basis of separation of multiple linked gene effects on mapping quantitative trait loci. Proc. Natl. Acad. Sci. USA 90: 10972–10976.[Abstract/Free Full Text]

ZENG, Z-B., 1994 Precision mapping of quantitative trait loci. Genetics 136: 1457–1468.[Abstract]

This article has been cited by other articles:

				L. Zhang, H. Li, Z. Li, and J. Wang Interactions Between Markers Can Be Caused by the Dominance Effect of Quantitative Trait Loci Genetics, October 1, 2008; 180(2): 1177 - 1190. [Abstract] [Full Text] [PDF]

				A. O. Bergland, A. Genissel, S. V. Nuzhdin, and M. Tatar Quantitative Trait Loci Affecting Phenotypic Plasticity and the Allometric Relationship of Ovariole Number and Thorax Length in Drosophila melanogaster Genetics, September 1, 2008; 180(1): 567 - 582. [Abstract] [Full Text] [PDF]

				R. D. Ball Quantifying Evidence for Candidate Gene Polymorphisms: Bayesian Analysis Combining Sequence-Specific and Quantitative Trait Loci Colocation Information Genetics, December 1, 2007; 177(4): 2399 - 2416. [Abstract] [Full Text] [PDF]

				S. Isobe, A. Nakaya, and S. Tabata Genotype Matrix Mapping: Searching for Quantitative Trait Loci Interactions in Genetic Variation in Complex Traits DNA Res, November 13, 2007; (2007) dsm020v1. [Abstract] [Full Text] [PDF]

				L. A. Robertson-Hoyt, C. E. Kleinschmidt, D. G. White, G. A. Payne, C. M. Maragos, and J. B. Holland Relationships of Resistance to Fusarium Ear Rot and Fumonisin Contamination with Agronomic Performance of Maize Crop Sci., September 1, 2007; 47(5): 1770 - 1778. [Abstract] [Full Text] [PDF]

				N. Yi, D. Shriner, S. Banerjee, T. Mehta, D. Pomp, and B. S. Yandell An Efficient Bayesian Model Selection Approach for Interacting Quantitative Trait Loci Models With Many Effects Genetics, July 1, 2007; 176(3): 1865 - 1877. [Abstract] [Full Text] [PDF]

				M. Zak, A. Baierl, M. Bogdan, and A. Futschik Locating Multiple Interacting Quantitative Trait Loci Using Rank-Based Model Selection Genetics, July 1, 2007; 176(3): 1845 - 1854. [Abstract] [Full Text] [PDF]

				N. Yi, S. Banerjee, D. Pomp, and B. S. Yandell Bayesian Mapping of Genomewide Interacting Quantitative Trait Loci for Ordinal Traits Genetics, July 1, 2007; 176(3): 1855 - 1864. [Abstract] [Full Text] [PDF]

				J. Szyda and J. Komisarek Statistical Modeling of Candidate Gene Effects on Milk Production Traits in Dairy Cattle J Dairy Sci, June 1, 2007; 90(6): 2971 - 2979. [Abstract] [Full Text] [PDF]

				B. Stich, J. Yu, A. E. Melchinger, H.-P. Piepho, H. F. Utz, H. P. Maurer, and E. S. Buckler Power to Detect Higher-Order Epistatic Interactions in a Metabolic Pathway Using a New Mapping Strategy Genetics, May 1, 2007; 176(1): 563 - 570. [Abstract] [Full Text] [PDF]

				H. Li, G. Ye, and J. Wang A Modified Algorithm for the Improvement of Composite Interval Mapping Genetics, January 1, 2007; 175(1): 361 - 374. [Abstract] [Full Text] [PDF]

				S. Viitala, J. Szyda, S. Blott, N. Schulman, M. Lidauer, A. Maki-Tanila, M. Georges, and J. Vilkki The Role of the Bovine Growth Hormone Receptor and Prolactin Receptor Genes in Milk, Fat and Protein Production in Finnish Ayrshire Dairy Cattle Genetics, August 1, 2006; 173(4): 2151 - 2164. [Abstract] [Full Text] [PDF]

				A. Baierl, M. Bogdan, F. Frommlet, and A. Futschik On Locating Multiple Interacting Quantitative Trait Loci in Intercross Designs Genetics, July 1, 2006; 173(3): 1693 - 1703. [Abstract] [Full Text] [PDF]

M. Perez-Enciso Multiple association analysis via simulated annealing (MASSA) Bioinformatics,