## Abstract

Many statistical methods have been developed to map multiple quantitative trait loci (QTL) in experimental cross populations. Among these methods, multiple-interval mapping (MIM) can map QTL with epistasis simultaneously. However, the previous implementation of MIM is for continuously distributed traits. In this study we extend MIM to ordinal traits on the basis of a threshold model. The method inherits the properties and advantages of MIM and can fit a model of multiple QTL effects and epistasis on the underlying liability score. We study a number of statistical issues associated with the method, such as the efficiency and stability of maximization and model selection. We also use computer simulation to study the performance of the method and compare it to other alternative approaches. The method has been implemented in QTL Cartographer to facilitate its general usage for QTL mapping data analysis on binary and ordinal traits.

MAPPING quantitative trait loci (QTL) is important for studying the genetic basis of quantitative trait variation. A number of statistical methods have been developed over the years for QTL mapping data analysis in designed experiments, such as those of Lander and Botstein (1989), Haley and Knott (1992), Jansen (1993), Zeng (1993, 1994), Sillanpää and Arjas (1998), and Kao *et al.* (1999). However, many of these statistical methods focus on continuous data. Ordinal traits are also common in many QTL mapping studies. These traits take values in one of several ordered categories. In quantitative genetics, we usually use a threshold model to model the genetic basis of binary and ordinal traits (Wright 1934a,b; Falconer 1965; Falconer and Mackay 1996). In this model, we assume that the categorical observation of a binary or ordinal trait is a reflection of an underlying continuously distributed liability subject to a series of thresholds that categorize phenotypes. Effects of QTL on observed phenotypes are modeled through the liability.

A number of studies have used this threshold model for QTL mapping analysis on binary and ordinal traits. Hackett and Weller (1995) and Xu and Atchley (1996) first studied a QTL mapping method for binary/ordinal traits based on composite interval mapping (CIM) (Zeng 1994). Visscher *et al.* (1996) compared the performance and statistical power of using a linear regression and a generalized linear model directly on a binary trait for QTL mapping analysis and observed that the two methods give quite similar results in detecting QTL and estimating QTL position. Yi and Xu (1999a,b) studied a few statistical issues for mapping QTL on binary traits in outbred populations. Broman (2003) proposed a method to deal with data with a spike in the trait distribution. In particular, Yi and Xu (2000, 2002) and Yi *et al.* (2004) reported a series of studies using a Bayesian approach for mapping QTL on binary and ordinal traits and studied several strategies for model selection in a Bayesian framework.

A major difference between an ordinal trait and a continuous trait is the number of trait values that a quantitative character may take: a few, say 2–10, for an ordinal trait (2 for a binary trait) and theoretically infinite for a continuous trait. As a result of this difference, it is more complicated to map QTL on an ordinal trait, since there is less information carried by the data. Therefore, it is important to use appropriate statistical methods that take the trait distribution into account for mapping QTL, particularly for mapping multiple QTL. For mapping multiple QTL, Kao *et al.* (1999) and Zeng *et al.* (1999) developed a method that fits a multiple-QTL model including epistasis on a trait and simultaneously searches the number, positions, and interaction of QTL. This method, called multiple-interval mapping (MIM), is based on maximum likelihood and combined with a model selection procedure and criterion. Compared with interval mapping (IM) (Lander and Botstein 1989) and CIM (Zeng 1994), MIM has a number of advantages, such as the improved statistical power in detecting multiple QTL (Zeng *et al*. 2000), facilitation for analyzing QTL epistasis, and coherent estimation of overall QTL parameters.

In this study, we extend MIM for mapping QTL on ordinal traits and study many associated statistical issues. The method is based on a threshold model, implemented in the framework of MIM and targeted to experimental populations, such as backcross and F_{2}. After introducing the models, we focus our discussion on many statistical issues, such as maximizing the likelihood function and model selection process. We also use simulations to investigate a few questions associated with analyzing multiple QTL on ordinal traits.

## METHODS

#### Threshold model and liability:

An imperative step in mapping QTL is to use appropriate models to connect trait values with QTL genotypes. For continuous data, MIM uses models that are described in *Genetic and statistical models* below. But for ordinal data, these models are not appropriate to be applied directly. However, with the help of a threshold model (Wright 1934a,b; Falconer 1965; Falconer and Mackay 1996), we can extend the models and methodology of MIM to ordinal data. The threshold model assumes that there is an underlying unobserved trait value, called liability, for the observed ordinal trait. The liability may be continuous. When it reaches a certain threshold, a categorical phenotype is observed. Thus, we can relate ordinal trait values to QTL genotypes by relating the ordinal data to their continuous liability first by the threshold model and then relating the liability to QTL genotypes by the regular genetic and statistical models.

Suppose in an experiment, *n* ordinal-scaled trait values are observed and are coded as 0, 1, … , *n* − 1. In addition, suppose *N* individuals are sampled for study. For the *i*th individual, let *z _{i}* be its ordinal-scaled trait value and

*y*its underlying liability, where

_{i}*i*= 1, … ,

*N*. By definition,

*z*takes a value from {0, 1, … ,

_{i}*n*− 1} and

*y*is from an unknown continuous distribution (the liability). These two values are related by the threshold model in the following way,where “⇔” represents “is equivalent to,”

_{i}*s*is a value from {0, 1, … ,

*n*− 1}, and γ

_{s}'s (

*s*= 0, 1, … ,

*n*− 1) are a set of fixed (unknown) values in an ascending order and are called

*thresholds*with γ

_{0}= −∞ and γ

_{n}= ∞. Briefly, the above relationship indicates that when the liability of an individual falls between γ

_{s}and γ

_{s+1}, its phenotypic value is

*s*; and on the other hand, if its phenotypic value is

*s*, its liability must fall between γ

_{s}and γ

_{s+1}.

#### Genetic and statistical models:

As mentioned earlier, mapping QTL requires a connection between phenotypic values and QTL genotypes. By using the threshold model, the observed phenotypic values are connected to the underlying (continuous) liability. The next step is to connect the liability with QTL genotypes. This can be done by using the usual genetic and statistical models that have been used in many previous studies. The statistical model is used to characterize the relationship between the liability of an ordinal trait and its components, which include a genotypic part determined by QTL genotype and a random variation part caused by environment. The genetic model is used to compute the genotypic value for an individual on the basis of its QTL genotype. Consider a trait determined by *m* diallelic QTL. On the basis of the partition of variance (Fisher 1918), the genetic model for a genotypic value *G* includes additive and dominant effects and interactions among loci. Specifically, the genotypic value for the *i*th individual can be expressed as Equation 1 for a backcross design and as Equation 2 for an F_{2} design (ignoring trigenic or higher-order interactions),(1)(2)where μ_{G} is the overall mean of genotypic values, *a _{j}* is the main effect of QTL

*j*in a backcross design or additive effect in an F

_{2}design, and

*d*is the dominant effect of QTL

_{j}*j*. In addition, (

*aa*)

_{jk}, (

*ad*)

_{jk}, and (

*dd*)

_{jk}are, respectively,

*additive*×

*additive*,

*additive*×

*dominant*, and

*dominant*×

*dominant*interaction effects between QTL

*j*and

*k. x*and

_{ij}*u*are the corresponding variables for the additive and dominant effects. With

_{ij}*Q*and

_{j}*q*representing alleles in the two inbred parental lines,

_{j}*x*takes values of for

_{ij}*Q*and for

_{j}Q_{j}*q*in a backcross design; and in an F

_{j}q_{j}_{2}design,and

With this specification of genetic model, the statistical model can be defined by(3)where *e _{i}* is usually assumed to be independently normally distributed with mean zero and variance σ

^{2}.

#### Likelihood analysis:

Given the genetic and statistic models, proper statistical methods are needed to obtain estimates for QTL parameters. One way is to find a set of parameter values that yield the highest probability for the observed data given the models. The maximum-likelihood (ML) method is designed for this purpose. To use the ML method, two steps are needed: deriving a likelihood function and then maximizing it using a reliable and efficient algorithm. The likelihood function is defined as the joint probability of the sample given the model. Maximum-likelihood analysis has been used in interval mapping (Lander and Botstein 1989), composite-interval mapping (Zeng 1993), and multiple-interval mapping (Kao *et al*. 1999). We describe the likelihood function for ordinal data in this section and show how to maximize it in the next section, using a backcross design as an example.

The likelihood function for ordinal data is defined as , where **Z** represents the phenotypes (in an ordinal scale), **M** the marker genotypes, the QTL position parameters (for example, measured as genetic distance from one end of a chromosome), the threshold model parameters γ_{s} (*s* = 1, … , *n* − 1), and the QTL effects (both the main and epistatic effects, such as *a*'s, *d*'s and *w*'s) on the underlying liability. In addition, let **I**(*S* = *s*) be a half-open interval bracketed by γ_{s} and γ_{s+1} as (γ_{s}, γ_{s+1}] (*s* = 0, … , *n* − 1) and **Q**_{ih} be the *h*th (*h* = 1, … , 2* ^{m}*) possible QTL genotype for the

*i*th individual. With the assumption of independent sampling, the likelihood function can be written as(4)where represents product,

**M**

_{i}is the marker genotype of the

*i*th individual, and

*P*(*|*)'s [such as ] are conditional probabilities that are explained below along with their relationship with the likelihood function.

is the probability for individual *i* having QTL genotype **Q**_{ih} given its marker genotype **M**_{i} and QTL positions . Formulas for computing this probability have been given in Kao and Zeng (1997) and for cases with missing marker genotypes, in Jiang and Zeng (1997).

is the probability for individual *i* having *y _{i}* for the underlying liability given its QTL genotypes and QTL effects .

is the probability of observing *z _{i}* given underlying liability

*y*and thresholds . It is one when

_{i}*y*is between and and zero otherwise. In other words, if , and if

_{i}*y*∉

_{i}**I**(

*S*=

*z*).

_{i}Define to be the cumulative distribution function (cdf) for . Note that . On the basis of Fubini's theorem and the properties of the cdf, we have

#### Parameter estimation:

Once a likelihood function is given, estimates of parameters can be made by finding a set of parameter values that maximize the likelihood function. Estimates obtained in this way are called maximum-likelihood estimates (MLEs). In our study, MLEs are obtained by maximizing a function, which is is the expected log-likelihood function of the complete data (Dempster *et al*. 1977), using a combined Newton-Raphson (NR)–EM algorithm (see appendixes a and b).

When QTL positions are selected (see more in the *Model selection* below), parameters such as QTL effects can be estimated using an approach combining NR and EM algorithms (see appendix b for the iterative process). This approach is investigated by X.-J. Qin and Z-B. Zeng (unpublished data). It is useful when the matrix **H** (the second derivative matrix) is not positive definite or close to nonpositive definite, under which the NR algorithm may break down (*i.e*., fail to converge). The approach is similar to the regular NR method except that a check point is added to examine whether the Cholesky decomposition succeeds. This consequently determines whether the course of the iterative process is kept on the NR algorithm or is redirected to the EM algorithm.

#### Model selection:

To search and select QTL, we adapt the MIM procedure of Kao *et al.* (1999) and Zeng *et al.* (1999). This is a stepwise model adaptation procedure combined with an initial model selection by markers (Zeng *et al*. 1999). The idea is to use a computationally more efficient procedure, such as stepwise marker selection, first to select an initial model and then to use several model modification procedures under the MIM model to optimize the model selection.

We use a stepwise logistic regression to select significant markers as an initial model (SAS Institute 1999). We recommend using a backward stepwise selection procedure with the significance level α = 0.01 or α = 0.05 for the *F*-statistic, if there are more samples than the number of markers; otherwise, a forward stepwise selection may be used.

After the initial model selection, the following procedure can be used to update the model:

Optimize QTL position estimation. The position of each QTL is updated by scanning the genomic region flanked by the two adjacent QTL to choose the best position as a new estimate, while fixing the other QTL at their current positions. This is performed for each QTL sequentially.

Search for a new QTL. The best position for a new putative QTL is first searched in the genome. The Bayesian information criterion (BIC) of the new model with one more QTL is compared with that of the previous model to decide whether to add this QTL in the model. If the QTL is added, the number of QTL is increased by one; otherwise the model is unchanged.

Test the current QTL effects. Each QTL effect is tested by comparing BICs of the models with or without the QTL effect conditional on other QTL. If some QTL are not significant, the number of QTL is reduced; otherwise the model is unchanged.

This procedure can be used iteratively until the model is unchanged. Usually the epistatic effects of QTL are searched and tested afterward among the QTL identified.

## SIMULATIONS AND RESULTS

We use computer simulations to investigate the performance of our approach. In each simulation (except indicated otherwise), 100 data sets are generated by Windows QTL Cartographer (Wang *et al*. 2005). Each data set includes 200 individuals and has one to eight chromosome(s). On each chromosome, 10 evenly distributed markers are simulated with 10 cM between adjacent markers. Various numbers of QTL are simulated, which are specified in Tables 1 and 8. For simplicity, all QTL have the same main effects. Backcross design is used for illustration. Each individual has a simulated liability score that is transformed, on the basis of the preset incidence rates, to a phenotype in binary/ordinal scale. For the purpose of comparison, all data sets are analyzed in three ways. We use the MIM module in QTL Cartographer to analyze the liability score (denoted by QTLC) for comparison. The binary/ordinal phenotype is analyzed either by our new method (denoted by bMIM) and by the MIM module in QTL Cartographer (denoted by QTLB), which ignores the fact that the phenotypic value is binary/ordinal. Different notations are used to represent different parameter setups to avoid potential confusion, such as 1C1Q for one-chromosome one QTL. Statistics are collected for different analysis methods on the percentage of data sets that obtain the correct number of QTL, the mean number of QTL detected, and the mean of QTL position estimates.

We use these simulations to discuss several questions, such as empirical critical values, effects of various factors on mapping results, suitability of using QTL Cartographer/MIM on binary traits directly, loss of information in mapping when data are scored in binary/ordinal scale as compared to a continuous distribution, epistasis, limitations of the new method, and estimation of heritability *h*^{2} for ordinal traits.

#### Empirical critical values:

Although the criteria and critical values used for model selection are very important issues in QTL mapping, they are not the main study subject in the current investigation. Therefore, no comprehensive investigation of model selection criteria is performed. Instead, we run a few sets of simulations to illustrate how critical values likely behave in mapping binary data for three situations: one QTL *vs.* two QTL, four QTL *vs.* five QTL, and eight QTL *vs.* nine QTL.

Critical values are estimated in two ways: direct data simulation (SM) and residual bootstrapping (RB) (Zeng *et al*. 1999). Both use a likelihood-ratio test statistic to test for the significance to add a QTL. For each case, 1000 test statistics are obtained under the null hypothesis. The 100(1 − α)th percentile of these statistics is chosen as the critical value at a significance level of α.

For SM, 1000 independent data sets are simulated on the basis of the given QTL positions and effects (equivalent to 1000 different experiments). For a specific test condition, say four QTL *vs.* five QTL, a four-QTL model is first established, by finding the set of four chromosome positions that yield the greatest likelihood among all combinations of four positions. A fifth position maximizing the likelihood given the four-QTL model is found. The test statistic is chosen as the likelihood ratio between the four-QTL model and the model with the fifth position.

For RB, all 1000 data sets are derived from one single data set, as outlined below. Again using four QTL *vs.* five QTL as an example, a data set is first simulated on the basis of the given QTL parameters. A four-QTL model is established as previously described. Denote as the estimated QTL positions and **b** as the estimated QTL effects. The *i*th individual then has a genotypic value of **X**_{ij}**b** with a probability of and has a phenotypic value 0 with a probability of . A new data set is generated by assigning the trait value of the *i*th individual to be 0 with a probability of *w*_{0} while keeping its marker genotype. The test statistic for the new data set is obtained using the same procedure as for SM.

Results from these two procedures are shown in Table 2. For few numbers of QTL and low values of *h*^{2} such as one QTL or two QTL with *h*^{2} = 0.1 or *h*^{2} = 0.3, two procedures obtain similar results. For higher values of *h*^{2} and greater numbers of QTL, more differences are seen in the results. For RB with the same heritability, similar results are obtained for testing conditions and , although they are quite different from results for . For different heritability, greater differences among results are seen. This trend has been seen before in cases for continuous traits by S. Wang and Z-B. Zeng (unpublished results). That is, heritability has some effects on critical values, and critical values change when the number of testing QTL is low and tend to be stable with increasing numbers of testing QTL. For SM, both testing conditions and heritability may greatly affect results. No clear trend is detected for SM. In addition, we explored an approach proposed by Lin and Zou (2004), which uses a score-type statistic and results in less variation among the empirical critical values (data not shown). Nevertheless, the inconsistency among critical values makes it difficult to use them in data analyses. Therefore, before a comprehensive investigation of this issue is complete, BIC will be used as a temporary solution in model selection, since this criterion has been used before and reasonable results are obtained (Zeng *et al*. 1999; Broman and Speed 2002).

#### Effects of various factors:

In real experiments, the range of parameter values varies dramatically. To understand what effects parameter value differences may have on mapping results is helpful in choosing appropriate mapping methods and in evaluating QTL mapping results. Two parameters that are likely to change from one experiment to another are investigated. They are the proportion of each category and heritability. Proportions of different categories may affect mapping results. For example, when categories are divided unevenly, few individuals may exist within a specific category that carries inadequate information for mapping QTL. With other conditions being the same, a greater heritability, which measures the proportion of total variation explained by genetic factors, usually means relatively larger QTL effects. Therefore, it should be easier in finding QTL when heritability increases.

Effects of proportions of categories are studied for binary cases with 20, 30, and 50% of individuals having phenotypic values of 0 (others having 1). Simulation is done under 4C4Q with *h*^{2} = 0.3. Numbers of data sets detecting various numbers of QTL, mean numbers of detected QTL, and estimated *h*^{2} (with its standard deviation) are shown in Table 3. Mapping results are also obtained for the same data sets using QTLB. For moderately unevenly divided data (such as a data set with 20% of individuals having a trait value of zero), the majority of data sets detect 2 or 3 QTL (totally 66%). For data sets divided more evenly, more data sets detect 2, 3, or 4 QTL (>90%). The mean number of detected QTL has increased from 2.57 for data with 20% zeros to ∼3 for cases with more evenly divided categories. The differences may be partially because with varying proportions of categories, the chances for individuals with various genotypes being in the same category are also changed.

Effects of heritability can be seen in several tables, such as Tables 4–7⇓⇓⇓. As expected, when *h*^{2} is higher, all methods obtain better mapping results given other situations being the same.

#### Suitability of using QTL Cartographer/MIM on ordinal data directly:

When the number of categories for an ordinal trait is relatively large, the data can be analyzed by approaches implemented for continuous traits, such as in Visscher *et al.* (1996), binary traits are analyzed directly by using a linear regression method proposed by Haley and Knott (1992). However, the use of Haley–Knott approximation can yield a substantially inflated residual variance (Xu 1995) and potentially decrease the power of detecting QTL. Here we use the maximum-likelihood-based method (MIM module in QTL Cartographer) to analyze binary/ordinal traits directly, with the hope of remedying part of the inflation.

Simulations are performed for different combinations of QTL and chromosomes (1C1Q, 2C2Q, 4C4Q, and 8C8Q) with various values of heritability (0.1, 0.3, 0.5, and 0.8). The results are shown in Tables 4–7⇑⇑⇑, with likelihood-ratio profiles from different approaches for 8C8Q shown in Figure 1. By comparing results between QTLB and QTLC and between QTLB and bMIM, we investigated whether/when QTLB can be used. Relative to QTLC, the efficiency of detecting QTL by QTLB increases with higher *h*^{2} and a lower number of QTL: from ∼60% for 4C4Q (8C8Q) with *h*^{2} = 0.3, under which the mean number of detected QTL is 1.87 (1.82) for QTLB and 2.87 (2.77) for QTLC, to almost 100% for 1C1Q with *h*^{2} = 0.3. Compared to bMIM, QTLB yields similar results when *h*^{2} is high and the number of QTL is low, but has a lower efficiency for low *h*^{2} or low *h*^{2} with a high number of QTL. For example, for 2C2Q with *h*^{2} = 0.1, ∼30% of data sets detect two QTL when bMIM is applied and only ∼10% when QTLB is used. Generally speaking, QTLB yield reasonable mapping results when *h*^{2} is high (say >0.4) and the number of QTL is low (say <4).

#### Loss of information:

Binary/ordinal data carry less information than continuous data due to at least two factors. One is that phenotypic values cannot be ranked in detail for ordinal data with only several categories. This lowers the resolution of QTL mapping and reduces the ability of finding QTL with small effects. The other is related to the shape of the distribution of phenotypic values. For ordinal data, trait values concentrate on several separate points instead of covering a region for a continuous trait. This may limit the ability to evaluate mapping results in terms of power and error rate. Since some traits may yield binary/ordinal data only for technical and practical reasons, investigating the efficiency of mapping QTL using binary/ordinal data relative to using continuous data can help us to better understand the limitation of analyzing QTL in experiments producing binary/ordinal data.

Comparing results from bMIM and QTLC in Tables 4–7⇑⇑⇑, we find that the efficiency of bMIM increases with higher *h*^{2}. For example, for 2C2Q with *h*^{2} = 0.1, 51 and 65% of position estimates by bMIM fall within 15 cM of the two corresponding simulated QTL, respectively, relative to 69 and 88% by QTLC; when *h*^{2} increases to 0.5, the results are 98% and 98% for bMIM and 100% and 100% for QTLC. However, the estimates for QTL numbers are relatively closer: for *h*^{2} = 0.1, 0.3, and 0.5, the numbers are 1.06, 1.93, and 2.16 for bMIM, compared to 0.90, 1.96, and 2.05 for QTL, respectively. A similar trend is also seen in multiple QTL and multiple chromosomes cases.

#### Epistasis:

To study epistasis, we adapt an approach using a stepwise model selection scheme in MIM, as described in Kao *et al.* (1999) and Zeng *et al.* (1999). In this scheme, two types of epistasis can be searched. The first type is the interaction between QTL with main effects. The second one is between QTL with main effects and QTL without main effects. For either type of interaction, epistatic effects are tested for statistical significance and are added to or dropped from the model on the basis of the testing results. Two case studies are used to illustrate and test our implementation. Details of parameter values are listed in Table 8. Briefly, for all cases three QTL are simulated on two chromosomes at the following positions: 25 cM (QTL 1) and 75 cM (QTL 2) on chromosome 1 and 35 cM (QTL 3) on chromosome 2. The interaction is between QTL pair 1–3 for all cases. For case 1, all QTL have main effects with *V*_{i}/*V*_{a} = 0.1 (case I-1) or 0.3 (case I-2), where *V*_{i} and *V*_{a} are the epistatic and additive variances, respectively. For case 2 (case II), QTL 1 and 2 have main effects, QTL 3 does not (hence QTL 3 enters into the model only through interaction), and *V*_{i}/*V*_{a} = 0.3.

Results are shown in Figure 2, based on 1000 simulated data sets with 300 individuals for each set of parameter values. Test statistics of the interaction among different QTL pairs are shown in Figure 2, a (for case I-1) and b (for case I-2). QTL pair with real interaction (QTL pair 1–3) is detected with a higher probability than other QTL pairs (QTL pairs 1–2 and 2–3). Indeed for 98% of case I-1 and for 88% of case I-2, QTL pair 1–3 has the largest statistic. In Figure 2c, the distribution of the case II test statistic, which is for interaction between detected QTL and regions without detected QTL, is shown. The distribution has a mean of ∼15 and a variance of ∼37. About 97% of the test statistics are significant when chi-square distribution is used as an approximation under the null hypothesis. In addition, 97% of these significant statistics have their corresponding QTL being in accordance with simulated data. Therefore, 94% of tests recover the simulated interaction and consequently, the average number of QTL detected changes from 2.08 when no epistasis is considered to ∼3.02 when epistasis is considered. The counts of detected QTL based on their estimated locations are shown in Figure 2d.

#### Limit:

It is also interesting to study parameters such as the minimal effect of QTL or maximal number of QTL that can be detected. This is useful in determining the applicability of a method and the validity of its results. Computer simulations are again used for the investigation. A preliminary result for a range of *h*^{2}-values and the number of QTL are shown in Table 9. It can be seen that the mean number of detected QTL increases when *h*^{2} increases. For example, for bMIM, the number increases from 0.17 to 0.61 for 1C1Q and from 0.18 to 0.37 for 2C2Q, respectively, when *h*^{2} changes from 0.01 to 0.05. This is partially expected because QTL effects are smaller when *h*^{2} is smaller and the number of QTL stays the same. However, no trend for the mean number of the detected QTL is seen when the same value of *h*^{2} and different numbers of simulated QTL are considered: the results fluctuate from 0.61 to 0.37 to 0.63 for 1C1Q, 2C2Q, and 4C4Q with *h*^{2} = 0.05, when bMIM is used. Another useful measurement is the percentage of detected QTL to simulated QTL (or the ratio between the mean number of detected QTL and the number of simulated QTL). For 1C1Q with *h*^{2} = 0.03, ∼30% of QTL could be detected by all three approaches; for four QTL with *h*^{2} = 0.1, percentages are <30. Note that the numbers are higher for bMIM for 4C4Q and 8C8Q. This may be due to lower critical values used in bMIM than they should be. Generally speaking, we expect that when QTL effects are <0.10, the percentage of detected QTL will be around or <5%, which is close to the rate of random errors and therefore may be close to QTL detection limits for these methods.

#### Approximation of *h*^{2}:

*R*^{2} of the fitted models can be used to approximate *h*^{2}, such as in QTLC and QTLB. For ordinal data using bMIM, the estimate of *h*^{2} can be approximated by an alternative form of *R*^{2} designed for logistic regression, suggested by Nagelkerke (1991),which is an adjusted form of proposed by Maddala (1983), Cox and Snell (1989), and Magee (1990),where *L*_{0} is the likelihood under the null model (no QTL), *L*_{1} is the maximum likelihood under the alternative model (a certain number of QTL exist), and *N* is the sample size. is the maximal possible value for and is equal to 1 − . Note that ranges between 0 and 1.

Using the simulated data sets of Tables 4–7⇑⇑, results of approximating *h*^{2} are summarized in Table 10. These results suggest that better approximation of *h*^{2} is obtained when underlying heritability (denoted by ) increases for a specific combination of numbers of QTL and chromosomes and that for the same value, a smaller number of QTL generally result in a better approximation of *h*^{2}. For example, for , estimates of *h*^{2} are 0.247, 0.270, and 0.301 by QTLB, bMIM, and QTLC, respectively, for 1C1Q, and 0.132, 0.198, and 0.205, respectively, for 4C4Q. These results are somehow expected: with a lower underlying heritability and/or greater number of QTL, QTL effects are smaller, and this will make it more difficult to detect these QTL. The model then will explain a smaller amount of total variation. In addition, QTLC has the best estimates and bMIM yields better results than QTLB does, especially for large number of QTL/chromosome situations.

## IMPLEMENTATION IN QTL CARTOGRAPHER

We have implemented the procedures of categorical trait MIM (CT–MIM) in version 2.5 of Windows QTL Cartographer (Wang *et al*. 2005). Categorical traits are coded in integer value and can be input into the program in the same way as continuous traits. One can use the categorical trait analysis to perform interval mapping (CT–IM) and multiple-interval mapping (CT–MIM) analysis. One can also use the regular IM, CIM, and MIM, implemented for continuous traits, to analyze the data for comparison.

For model selection in CT–MIM, we implemented several interactive procedures. They include a procedure that selects an initial model on the basis of forward or backward stepwise logistic regression on markers, a procedure to search for more QTL one at a time, a procedure to search for epistasis of QTL, a procedure to test effects of selected QTL, a procedure to optimize positions of QTL, and a procedure to output the complete information of the selected model. These procedures can be used interactively in practical data analysis to explore the data and compare different models. For more information see the software manual (http://statgen.ncsu.edu/qtlcart/WQTLCart.htm).

## DISCUSSION

Kao *et al.* (1999) developed MIM for QTL analysis that fits a multiple-QTL model and simultaneously searches for the positions and interaction patterns of multiple QTL. This multiple-QTL-oriented approach has a number of advantages, such as improved statistical power in QTL detection, facilitation for analyzing QTL epistasis, and coherent estimation of overall QTL parameters. The implementation of the MIM method in QTL Cartographer (Basten *et al*. 2005) and Windows QTL Cartographer (Wang *et al*. 2005), both freely available at http://statgen.ncsu.edu/qtlcart/, has greatly facilitated applications of the method for general QTL mapping data analysis. However, the MIM method of Kao *et al.* (1999) and its implementation in QTL Cartographer and Windows QTL Cartographer is for continuous traits. Although extensive research has been made that takes multiple QTL into account in mapping analysis on binary and ordinal traits (*e.g.*, Yi and Xu 2000; Yi *et al*. 2004), no computer program is available for QTL mapping data analysis on ordinal traits under the MIM framework.

In this study, we extend MIM to ordinal traits on the basis of a threshold model. This extension utilizes the properties and advantages of MIM for QTL mapping analysis on ordinal traits. The method fits a model of multiple-QTL effects and epistasis on the underlying liability and searches for the number and positions of QTL and epistasis simultaneously. It has similar advantages to MIM on the statistical power of QTL detection and estimation of overall QTL parameters. The implementation of this method in Windows QTL Cartographer will greatly facilitate the general usage of the methodology for QTL mapping data analysis on binary and ordinal traits.

Using simulations, we investigated several statistical issues, such as the effect of trait distribution on QTL mapping results, comparison of QTL mapping on an ordinal trait and on a continuous trait, and the statistical power of the method for QTL detection. As expected, the larger the number of trait categories, the higher the statistical power for QTL detection. There is not much difference in mapping results to regard an ordinal trait with five or more categories as a continuous trait in QTL analysis (data not shown). Also it is interesting to observe that if we regard a binary trait as a continuous trait using the current MIM in QTL Cartographer, the mapping result is actually quite comparable to that using the threshold MIM model if the heritability is reasonably high. Of course, the threshold MIM model is always more powerful and appropriate for QTL mapping analysis on binary and ordinal traits.

In studying ordinal traits, the trait value of an individual may be misspecified due to measurement error. For a binary trait, Rousseeuw and Christmann (2003) used a “hidden logistic regression model” with the assumption that an observed response has a small chance to be measured with error. This can occur when a binary or ordinal phenotype is difficult to classify. A model with measurement error similar to that of Rousseeuw and Christmann (2003) can also be used for our analysis. This can be done by reassigning the value of in Equation 4. Namely, instead of being either 1 or 0, it can be 1 − ε_{i} and ε_{i}, where ε_{i} is a small nonnegative value for error rate. This error rate can be assumed to be the same for all observations or different for different observations.

In this article, we used the maximum-likelihood approach for mapping multiple QTL on binary and ordinal traits. The Bayesian approach has also been used extensively for QTL mapping analysis in designed experiments, such as in Thomas and Cortessis (1992), Hoeschele and van Raden (1993a,b), Satagopan and Yandell (1996), and Sillanpää and Arjas (1998, 1999). For binary and ordinal traits, a series of studies have been performed by Xu and Atchley (1996), Yi and Xu (1999a,b, 2000), and Yi *et al.* (2004) to develop statistical methods for mapping multiple QTL under a Bayesian framework combined with Markov chain Monte Carlo sampling with a reversible-jump algorithm for model selection. These methods, based on a threshold model for binary/ordinal traits, are comparable to our maximum-likelihood method. However, despite the extensive studies performed, no user-friendly software is publicly available for QTL mapping data analysis on binary/ordinal traits. The statistical methods described in this article will be implemented in QTL Cartographer and Windows QTL Cartographer and publicly distributed at http://statgen.ncsu.edu/qtlcart/ for general usage of mapping multiple QTL on binary and ordinal traits.

There are still some issues that deserve further investigation. Currently, we use a procedure suggested by Lin and Zou (2004) to estimate the threshold at each step of searching for new QTL to aid in model selection. This procedure performs a function similar to a permutation test, but is numerically much more efficient. However, it is still not quite clear yet what significance level one needs to use in this stepwise procedure in the context of model selection for multiple QTL with epistasis. We will further pursue this line of research.

## APPENDIX A: THE 𝒬-FUNCTION AND ITS DERIVATIVES

The -function is defined as the expectation of log-likelihood function for the complete data (Dempster *et al*. 1977). In our study, QTL genotypes are unknown and have to be inferred from marker genotypes. This can be dealt with by the -function and the process is briefly described below aswhere is the vector of parameters to be estimated, and the superscript *t* indicates the *t*th cycle of the iteration. Since the real **B** is unknown, we compute the -function on the basis of its value at the *t*th stage. This is symbolized as . The subscript for the expectation sign *E* indicates that the expectation is computed using a specific set of values. In addition, with missing data, the complete likelihood *L*_{c} has to be computed on the basis of observed data {**Z**, **M**}_{obs}. Using definitions of expectation and conditional probability and assuming independent sampling, we have(A1)where *C* is the probability of the observed data and is a constant. Therefore, *C* can be omitted where maximization of the -function is concerned. By Bayes' theorem and noting that *P*_{B}_{=B}^{(t)}(*z _{i}*|

**Q**

_{ih},

**M**

_{i})=

*P*

_{B}_{=B}

^{(t)}(

*z*|

_{i}**Q**

_{ih}), we have(A2)

Here, can be considered as a posterior probability for **Q**_{ih}. Equation A1 can then be written as(A3)

With a logistic-distributed liability and denoting , Equation A3 can be written as with(A4)where subscript *ih* indicates the *h*th QTL genotype for the *i*th individual, and is a parameter vector for the thresholds and QTL effects. In addition, , where **l**_{k} (*k* = 0, … , *n* − 1) is an *n* × 1 vector with all elements being 0 except the (*k* + 1)th element being 1. Defineand

Noting that ∂π_{k,ih}/∂**B**_{0} = π_{k,ih}(1 − π_{k,ih})**x**_{k,ih}, we have(A5)

Further differentiating the first derivative, we have the second derivative as(A6)where

## APPENDIX B: THE ITERATIVE PROCESS FOR THE NR–EM ALGORITHM

For a set of fixed QTL positions, an approach combining Newton–Raphson (NR) and EM algorithms can be used to obtain estimates of other parameters such as QTL effects. A brief description of the iterative process is given below:

Initialize .

Update in Equation A2 (the E-step) and π

_{ih}in Equation A4.Obtain the first derivative (

**g**^{(t)}) and second derivative (**H**^{(t)}) using Equations A5 and A6.Update using a formula

where the superscript “−1” indicates the inverse of a matrix, and the inverse of **H** could be obtained by Cholesky decomposition. Several EM steps may be performed in the case that the decomposition fails.

Find the new value of the -function at by Equation A3.

Determine whether the iteration process converges, usually by comparing the relative change of the to a convergence tolerance (δ). If the change of the -function is smaller than δ, stop; if not, go back to step 2.

Note that two parameters (β and δ) in the above process need to be preset. β is a scalar and characterizes the step length in the gradient direction. It should not be too large (missing maximum) or too small (slow convergence). During computation, β will start at one and reduce its value gradually in each cycle if needed, but it will be reset to one in a new iteration cycle. The value of δ can be determined through several methods. Here, we take the simplest one: set δ to a fixed small number, say 10^{−8}.

## Acknowledgments

We thank Yongqiang Tang for help in computer programming. This work was partially supported by National Institutes of Health grant GM45344 and U.S. Department of Agriculture Plant Genome grant 2003-00754.

## Footnotes

↵1

*Present address:*Department of Molecular Biology and Genetics, Cornell University, Ithaca, NY 14853.Communicating editor: M. S. McPeek

- Received December 12, 2005.
- Accepted March 22, 2006.

- Copyright © 2006 by the Genetics Society of America