### Multivariate linear model

Consider *n* individuals derived from a backcross population crossed from two inbred lines with observations on some densely distributed codominant markers and on *m* quantitative traits. Supposed that the maximum number of QTL is *p*, the phenotypic value *y*_{
ki
}of individual *i* for *k* th trait can be described by the following multivariate linear model:

{y}_{ki}={b}_{k0}+{\displaystyle \sum _{j=1}^{p}{\gamma}_{kj}{x}_{kij}{b}_{kj}}+{e}_{ki},

(1)

for *i* = 1, 2, ..., *n* and *k* = 1, 2, ..., *m*, where *γ*_{
kj
}is model indicator variable, indicating the *j* th QTL of *k* th trait included (1) or excluded (0) from the model; *b*_{k 0}is population mean; *b*_{
kj
}is QTL effect; *x*_{
kij
}is QTL genotype, if QTL genotype is homozygote *x*_{
kij
}= 1, otherwise -1; *e*_{
ki
}is residual error and assumed to follow multivariate normal distribution. If we denote equation (1) by matrix, it can be expressed as:

{y}_{i}={b}_{0}+{\displaystyle \sum _{j=1}^{p}{\mathbf{\Phi}}_{j}{X}_{ij}{b}_{j}}+{e}_{i},

(2)

for *i* = 1, 2, ..., *n*, where **y**_{
i
}= [*y*_{1i}, *y*_{2i}, ..., *y*_{
mi
}]^{T}, **b**_{0} = [*b*_{10}, *b*_{20}, ..., *b*_{m 0}]^{T}, **b**_{
j
}= [*b*_{1j}, *b*_{2j}, ..., *b*_{
mj
}]^{T}, **e**_{
i
}= [*e*_{1i}, *e*_{2i}, ..., *e*_{
mi
}]^{T}. They are all (1 × *m*) column vectors. Equation (3) is QTL genotype matrix and Equation (4) is model indicator matrix, they are all (*m* × *m*) diagonal matrix.

{X}_{ij}=\left[\begin{array}{cccc}{x}_{1ij}& 0& \cdots & 0\\ 0& {x}_{2ij}& \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0& 0& \cdots & {x}_{mij}\end{array}\right]

(3)

{\mathbf{\Phi}}_{j}=\left[\begin{array}{cccc}{\gamma}_{1j}& 0& \cdots & 0\\ 0& {\gamma}_{2j}& \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0& 0& \cdots & {\gamma}_{mj}\end{array}\right]

(4)

### Prior specification

The prior distribution of each QTL effect vector **b**_{
j
}is multivariate normal distribution, *p*(**b**_{
j
}) ~ *N*(0, {\mathbf{\Sigma}}_{{B}_{j}}), where {\mathbf{\Sigma}}_{{B}_{j}} is the hyper-parameter, and We take {\mathbf{\Sigma}}_{{B}_{j}}={\left[X{.}_{j}^{T}{\mathbf{\Sigma}}_{e}^{-1}X{.}_{j}\right]}^{-1}\cdot n, which is simply an extension from Bayesian single trait analysis [15]. The importance of the choice of the hyper-parameter will be discussed later. In a large backcross population and under the definition of *x*_{
mij
}(-1 or 1), {\mathbf{\Sigma}}_{{B}_{j}} can be simplified as {\mathbf{\Sigma}}_{{B}_{j}} = **Σ**_{
e
}. The prior of the covariance matrix of residual error follows Inverse Wishart distribution, **Σ**_{
e
}~ *Wishart*^{-1}(*v*_{
e
}, {S}_{e}^{2}), where, *v*_{
e
}and {S}_{e}^{2} are prior degree of freedom and covariance matrix of residual error, respectively, and can be obtained from other method, such as CIM based multitrait analysis [2], *etc*. The prior distribution of population mean **b**_{0} is normal distribution with mean and variance equal to those calculated by phenotypic values. The prior probability distribution of QTL position *λ*_{
kj
}is uniform distribution with bounds of two flanking markers, *p*(*λ*_{
kj
}) = 1/*d*_{
j
}, where *d*_{
j
}is length of the interval where *j* th QTL is confined. Assuming that epistatic effect is absent, the prior inclusion probability for *j* th effect can be expressed as *p*(*γ*_{
kj
}= 1) = 1 - *l*_{
k
}/*L*_{
k
}]^{1/N}(see also [15]), where *l*_{
k
}is the prior expected number of main-effect QTL, and could be roughly estimated with the use of standard genome scans; *N* is the number of possible main effects for each QTL and equal to 1 in BC family [15]; *L*_{
k
}is the upper bound of QTL number, and equals to the number of marker interval in our simulation study, while in another approach suggested by Yi [15]*L*_{
k
}is taken as 3 + 3·{\sqrt{l}}_{k}, which causes the model space to reduce dramatically [15].

### Joint posterior density

The observable variables include phenotypic values, y={\left\{{y}_{i}\right\}}_{i=1}^{n} and marker information, m={\left\{{m}_{ij}\right\}}_{i=1,j=1}^{n,p}. The unobservable variables include population mean, {b}_{0}={\left\{{b}_{k0}\right\}}_{k=1}^{m}; QTL effects, b={\left\{{b}_{j}\right\}}_{j=1}^{p}; QTL genotypes, X={\left\{{X}_{ij}\right\}}_{i=1,j=1}^{n,p}; model indicator variables, \mathbf{\Phi}={\left\{{\mathbf{\Phi}}_{j}\right\}}_{j=1}^{p}; (co)variance of residual error, **Σ**_{
e
}, and QTL positions, \mathbf{\lambda}={\left\{{\lambda}_{kj}\right\}}_{k=1,j=1}^{m,p}. Let **θ** be the vector of hyper-parameters, **Θ** = {**b**_{0}, **b**, **Σ**_{
e
}, **λ**, **X**, **Φ**}, then the joint prior density of the unobservable variables is denoted by *p*(**Θ**|**θ**). The joint posterior probability of **Θ**, given the observable variables **y** and **m**, can be expressed as:

*p*(**Θ**|**y**, **m**) ∝ *p*(**Θ**|**θ**)·*p*(**y**, **m**|**Θ**), (2)

where, *p*(**y**, **m**|**Θ**) is the likelihood and can be written as:

*p*(**y**, **m**|**Θ**) = *p*(**y**|**Θ**)·*p*(**m**|**Θ**), (6)

where *p*(**y**|**Θ**) is multivariate normal density, and *p*(**m**|**Θ**) can be derived from a Markov model [14].

### MCMC sampling

MCMC algorithm generates samples from Markov chains which converge to the posterior distribution of parameters, without the constant of proportionality being calculated. From these posterior samples, summary statistic of the posterior distribution can be calculated. MCMC algorithm proceeds as follows:

a. Initialize all parameters with values in their legal domain.

b. Update the population mean **b**_{0}.

c. Update the QTL effects vectors {\left\{{b}_{j}\right\}}_{j=1}^{p}.

d. Update the variance-covariance matrix **Σ**_{
e
}of the residual error.

e. Update the QTL genotype indicator matrices {\left\{{X}_{ij}\right\}}_{i=1}^{n} and the QTL location vectors {\left\{{\mathbf{\lambda}}_{kj}\right\}}_{k=1}^{m} jointly, for *j* = 1, 2,..., *p*.

f. Update the model indicator variable matrices {\left\{{\mathbf{\Phi}}_{j}\right\}}_{j=1}^{p}.

The conditional posterior distribution of the population mean **b**_{0} is multivariate normal with mean

{\overline{b}}_{0}={\left[{\displaystyle \sum _{i=1}^{n}({\mathbf{\Sigma}}_{e}^{-1})}\right]}^{-1}{\displaystyle \sum _{i=1}^{n}{\mathbf{\Sigma}}_{e}^{-1}({y}_{i}}-{\displaystyle \sum _{j=1}^{p}{\mathbf{\Phi}}_{j}{X}_{ij}{b}_{j}}),

(7)

and variance-covariance matrix

{\mathbf{\Sigma}}_{{b}_{0}}={\left[{\displaystyle \sum _{i=1}^{n}({\mathbf{\Sigma}}_{e}^{-1})}\right]}^{-1}.

(8)

The conditional posterior distribution of the QTL effect **b**_{
j
}is sampled from multivariate normal distribution with mean

{\overline{b}}_{j}={\left[{\mathbf{\Sigma}}_{B}^{-1}+{\displaystyle \sum _{i=1}^{n}({X}_{ij}^{T}{\mathbf{\Phi}}_{j}^{T}{\mathbf{\Sigma}}_{e}^{-1}{\mathbf{\Phi}}_{j}X{}_{ij})}\right]}^{-1}{\displaystyle \sum _{i=1}^{n}{X}_{ij}^{T}{\mathbf{\Phi}}_{j}^{T}{\mathbf{\Sigma}}_{e}^{-1}({y}_{i}-}{\displaystyle \sum _{j\ne 1}^{p}{\mathbf{\Phi}}_{j}{X}_{ij}{b}_{j}}-{b}_{0}),

(9)

and variance-covariance matrix

{\mathbf{\Sigma}}_{{b}_{j}}={\left[{\mathbf{\Sigma}}_{B}^{-1}+{\displaystyle \sum _{i=1}^{n}({X}_{ij}^{T}{\mathbf{\Phi}}_{j}^{T}{\mathbf{\Sigma}}_{e}^{-1}{\mathbf{\Phi}}_{j}X{}_{ij})}\right]}^{-1}.

(10)

The posterior distribution of the residual error follows inverted Wishart distribution,

{\mathbf{\Sigma}}_{e}~Wishar{t}^{-1}(d{f}_{e}+{\nu}_{e},{\mathbf{\Omega}}^{T}\mathbf{\Omega}+{S}_{e}^{2}),

(11)

where \mathbf{\Omega}={y}_{i}-{\displaystyle \sum _{j=1}^{p}{\mathbf{\Phi}}_{j}{X}_{ij}{b}_{j}}-{b}_{0} and *df*_{
e
}= *n*.

In step e, the QTL locations and QTL genotype matrices are updated jointly. For locus *j*, we can firstly sample a new QTL position for each trait from their prior distribution (described later), then sample the QTL genotype matrices {\left\{{X}_{ij}\right\}}_{i=1}^{n} on the new position using equation (15), and finally, they are updated by the efficient Metropolis-Hastings algorithm [20, 21]. Because the sampling of **X**_{
ij
}is too complicate and we are going to firstly describe it. Due to the QTL genotype *x*_{
kij
}has two possible values (-1 or 1) in BC line, if *m* traits are investigated jointly, **X**_{
ij
}has 2^{m}kinds of possible formations, and the general pattern of **X**_{
ij
}can be written as:

{H}_{ij,{z}_{1}{z}_{2}\cdots {z}_{m}}=\left[\begin{array}{cccc}{x}_{1ij}={z}_{1}& 0& \cdots & 0\\ 0& {x}_{2ij}={z}_{2}& \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0& 0& \cdots & {x}_{mij}={z}_{m}\end{array}\right],

(12)

where, *z*_{1}, *z*_{2}, ..., *z*_{
m
}∈ {-1, 1}. For clarity, we omit the subscript *ij* from {H}_{ij,{z}_{1}{z}_{2}\cdots {z}_{m}} and present formulas {H}_{{z}_{1}{z}_{2}\cdots {z}_{m}} to denote the genotype matrix of *i* th individual and *j* th loci. Because the QTL genotypes *x*_{
kij
}of *i* th individual in the *j* th interval for all traits may be correlated, the joint prior probability of the genotype matrix **X**_{
ij
}can't be simply expressed by the following equation:

\begin{array}{c}p({X}_{ij}={H}_{{z}_{1}{z}_{2}\cdots {z}_{m}}|{\mathbf{\lambda}}_{j},{m}_{i,j},{m}_{i,j+1})=p({x}_{1ij}={z}_{1},{x}_{2ij}={z}_{2},\cdots ,{x}_{mij}={z}_{m}|{\mathbf{\lambda}}_{j},{m}_{i,j},{m}_{i,j+1})\\ ={\displaystyle \mathbf{\prod}_{k=1}^{m}p({x}_{kij}={z}_{k}|{m}_{i,j},{m}_{i,j+1})}\end{array}.

(13)

Instead, it can be derived from the Markov model (see Equation 14), assuming that the order of markers and QTL is *M*_{
j
}*Q*_{1}*Q*_{2} ... *Q*_{
m
}*M*_{j+1}(see Figure 7), where, *Q*_{1}, *Q*_{2}, ..., and *Q*_{
m
}denote the QTL respectively affecting trait 1, trait 2, ..., and trait *m* in *j* th marker interval. Indicator variables *x*_{1ij}, *x*_{2ij}, ..., and *x*_{
mij
}denote the genotypes of these QTL.

\begin{array}{l}p({X}_{ij}={H}_{{z}_{1}{z}_{2}\cdots {z}_{m}}|{m}_{i,j},{\mathbf{\lambda}}_{j},{m}_{i,j+1})=p({x}_{1ij}={z}_{1},{x}_{2ij}={z}_{2},\cdots ,{x}_{mij}={z}_{m}|{m}_{i,j},{\mathbf{\lambda}}_{j},{m}_{i,j+1})\\ \begin{array}{c}=p({x}_{1ij}={z}_{1}|{m}_{i,j},{\mathbf{\lambda}}_{1j},{m}_{i,j+1})\cdot p({x}_{2ij}={z}_{2}|{m}_{i,j},{\mathbf{\lambda}}_{2j},{x}_{1ij},{m}_{i,j+1})\\ \phantom{\rule{0.1em}{0ex}}\times \cdot \cdots \cdot p({x}_{mij}={z}_{m}|{m}_{i,j},{x}_{1ij},{x}_{2ij},\cdots ,{x}_{(m-1)ij},{\mathbf{\lambda}}_{mj},{m}_{i,j+1}),\end{array}\end{array}

(14)

If no segregation interference is considered, the joint prior probability can be factorized into equation (14), and each term in equation (14) can be derived from Haldane map function. Only the first term in equation (14) is conditional on two flanking markers; others are not only conditional on two flanking markers but also on the genotypes of all the QTL prior to the interested one. If double recombination is ignored [2], each term in equation (14) can be inferred only by the genotype of the left nearest loci (marker or QTL) and the right marker, then equation (14) can be simplified as:

\begin{array}{l}p({X}_{ij}={H}_{{z}_{1}{z}_{2}\cdots {z}_{m}}|{m}_{i,j},{\mathbf{\lambda}}_{j},{m}_{i,j+1})=p({x}_{1ij}={z}_{1},{x}_{2ij}={z}_{2},\cdots ,{x}_{mij}={z}_{m}|{m}_{i,j},{\mathbf{\lambda}}_{j},{m}_{i,j+1})\\ =p({x}_{1ij}={z}_{1}|{m}_{i,j},{\mathbf{\lambda}}_{1j},{m}_{i,j+1})\cdot p({x}_{2ij}={z}_{2}|{x}_{1ij},{\mathbf{\lambda}}_{2j},{m}_{i,j+1})\\ \phantom{\rule{0.1em}{0ex}}\times \cdot \cdots \cdot p({x}_{mij}={z}_{m}|{x}_{(m-1)ij},{\mathbf{\lambda}}_{mj},{m}_{i,j+1}),\end{array}

(15)

Each term in equation (15) can be easily inferred.

It is worth mentioning that we assume the sequence of markers and QTL is *M*_{
j
}*Q*_{1}*Q*_{2} ... *Q*_{
m
}*M*_{j+1}, and in fact, the sequence of QTL may be variable in each round of updating. Therefore, we should firstly ascertain the sequence in each round, and then construct the appropriate formula to calculate the joint prior probability of the QTL genotype *p*(**X**_{
ij
}= {H}_{{z}_{1}{z}_{2}\cdots {z}_{m}}|*m*_{i,j},**λ** *j*,*m*_{i,j+1}) according above rules. For clarity, we take an example to demonstrate it. Consider 3 QTL *Q*_{1}, *Q*_{2}, and *Q*_{3} that affect 3 traits respectively in an interval. Assuming that in a certain round the sequence of markers and QTL is *M*_{
j
}*Q*_{3}*Q*_{1}*Q*_{2}*M*_{j+1}, then the formula for calculating the joint prior probability of the QTL genotype can be written as:

\begin{array}{l}p({X}_{ij}={H}_{{z}_{1}{z}_{2}{z}_{3}}|{m}_{i,j},{\mathbf{\lambda}}_{j},{m}_{i,j+1})=p({x}_{1ij}={z}_{1},{x}_{2ij}={z}_{2},{x}_{3ij}={z}_{3}|{m}_{i,j},{\mathbf{\lambda}}_{j},{m}_{i,j+1})\\ \begin{array}{c}=p({x}_{3ij}={z}_{3}|{m}_{i,j},{\mathbf{\lambda}}_{3j},{m}_{i,j+1})\cdot p({x}_{1ij}={z}_{1}|{x}_{3ij},{\mathbf{\lambda}}_{1j},{m}_{i,j+1})\\ \times p({x}_{2ij}={z}_{2}|{x}_{1ij},{\mathbf{\lambda}}_{2j},{m}_{i,j+1}).\end{array}\end{array}

Once we obtain the joint prior probability of the QTL genotype, the joint conditional posterior probability of **X**_{
ij
}can be expressed as:

p({X}_{ij}={H}_{{z}_{1}{z}_{2}\cdots {z}_{m}}|{y}_{i},\cdots )=\frac{f({y}_{i}|{X}_{ij}={H}_{{z}_{1}{z}_{2}\cdots {z}_{m}},\cdots )p({X}_{ij}={H}_{{z}_{1}{z}_{2}\cdots {z}_{m}}|{\lambda}_{j},{m}_{ij},{m}_{i,j+1})}{{\displaystyle {\sum}_{{h}_{1}=-1}^{1}{\displaystyle {\sum}_{{h}_{2}=-1}^{1}\cdots {\displaystyle {\sum}_{{h}_{m}=-1}^{1}f({y}_{i}|{X}_{ij}={H}_{{h}_{1}{h}_{2}\cdots {h}_{m}},\cdots )p({X}_{ij}={H}_{{h}_{1}{h}_{2}\cdots {h}_{m}}|{\lambda}_{j},{m}_{ij},{m}_{i,j+1})}}}}

(16)

where f({y}_{i}|{X}_{ij}={H}_{{z}_{1}{z}_{2}\cdots {z}_{m}},\cdots ) is likelihood, and follows multivariable normal distribution,

f({y}_{i}|{X}_{ij}={H}_{{z}_{1}{z}_{2}\cdots {z}_{m}},\cdots )=\frac{1}{{(2\pi )}^{m/2}{\left|{\mathbf{\Sigma}}_{e}\right|}^{1/2}}\mathrm{exp}\phantom{\rule{0.1em}{0ex}}\left\{-\frac{1}{2}{({y}_{i}-{\displaystyle \sum _{j=1}^{p}{\mathbf{\Phi}}_{j}{X}_{ij}{b}_{j}}-{b}_{0})}^{T}{\mathbf{\Sigma}}_{e}^{-1}({y}_{i}-{\displaystyle \sum _{j=1}^{p}{\mathbf{\Phi}}_{j}{X}_{ij}{b}_{j}-{b}_{0}})\right\}

(17)

Once we have calculated 2^{m}possible posterior probabilities for the corresponding QTL genotype matrices, we are going to sample one genotype matrix according to their posterior probabilities. We firstly constructed the cumulative probability function *F*(*d*) by accumulating the 2^{m}probabilities in an arbitrary sequence for *d* = 1, 2, ..., 2^{m}and *F*(0) = 0, which is a discrete distribution; then sampled a random number from uniform distribution, *u* ~ *U*[0,1]; and compared *u* with *F*(*d*), if *F*(*d* - 1) <*u* ≤ *F*(*d*), then the *d* th genotype matrix is accepted.

The new sampled QTL genotype matrices {\left\{{X}_{ij}\right\}}_{i=1}^{n} are only the proposal value, which should be updated along with the proposal QTL position vector **λ**_{
j
}= [*λ*_{1j}, *λ*_{2j}, ..., *λ*_{
mj
}] by the Metropolis-Hastings algorithm [20, 21]. For each trait, the new proposal position is sampled around the existing one from uniform distributions, {\mathbf{\lambda}}_{kj}^{\ast} ~ [*λ*_{
kj
}- *δ*, *λ*_{
kj
}+ *δ*), where *δ* is tuning parameter, usually taking a value of 1 or 2 cM. The new position vector is denoted by {\mathbf{\lambda}}_{j}^{\ast}=\left[{\mathbf{\lambda}}_{1j}^{\ast},{\mathbf{\lambda}}_{2j}^{\ast},\cdots ,{\mathbf{\lambda}}_{mj}^{\ast}\right]; then the new QTL genotype matrix {X}_{ij}^{\ast} is sampled conditionally on the new position using equation (16); finally, the position vector {\mathbf{\lambda}}_{j}^{\ast} and genotype matrices {\left\{{X}_{ij}\right\}}_{i=1}^{n} are accepted jointly with probability equal to min(1,α), where

\mathit{\alpha}=\frac{{\displaystyle {\prod}_{i=1}^{n}p({y}_{i}|{X}_{ij}^{\ast},{\mathbf{\lambda}}_{j}^{\ast},\cdots )p({X}_{ij}^{\ast}|{\mathbf{\lambda}}_{j}^{\ast},\cdots )p({\mathbf{\lambda}}_{j}^{\ast})}}{{\displaystyle {\prod}_{i=1}^{n}p({y}_{i}|{X}_{ij},{\mathbf{\lambda}}_{j},\cdots )p({X}_{ij}|{\mathbf{\lambda}}_{j},\cdots )p({\mathbf{\lambda}}_{j})}}\cdot \frac{q({X}_{ij}|{y}_{i},\cdots )q({\mathbf{\lambda}}_{j})}{q({X}_{ij}^{\ast}|{y}_{i},\cdots )q({\mathbf{\lambda}}_{j}^{\ast})},

(18)

*p*({\mathbf{\lambda}}_{j}^{\ast}) and *p*(**λ**_{
j
}) is the prior probability of new and old position respectively, and they are cancelled out under uniform prior distribution; p({X}_{ij}^{\ast}|{\mathbf{\lambda}}_{j}^{\ast},\cdots ) and *p*(**X**_{
ij
}|**λ**_{
j
}, ...) is the prior probability of QTL genotype conditional on new and old position, which has been described detailed previously; \frac{q({X}_{ij}|{y}_{i},\cdots )}{q({X}_{ij}^{\ast}|{y}_{i},\cdots )}=\frac{p({X}_{ij}|{y}_{i},\cdots )}{p({X}_{ij}^{\ast}|{y}_{i},\cdots )} and \frac{q({\mathbf{\lambda}}_{j})}{q({\mathbf{\lambda}}_{j}^{\ast})}=\frac{{\displaystyle {\prod}_{k=1}^{m}p({\mathbf{\lambda}}_{kj})}}{{\displaystyle {\prod}_{k=1}^{m}p({\mathbf{\lambda}}_{kj}^{\ast})}}, are all proposal ratio.

In step f, block sampling of the indicator variable matrix **Φ**_{
j
}is expected to have a better performance than separately updating each *γ*_{
kj
}in **Φ**_{
j
}. Due to there are two possible values (0 or 1) for each model indicator *γ*_{
kj
}, if *m* traits are investigated jointly, each model indicator matrix **Φ**_{
j
}has 2^{m}kinds of formations. The general formula of it can be written as:

{W}_{j,{w}_{1}{w}_{2}\cdots {w}_{m}}=\left[\begin{array}{cccc}{\gamma}_{1j}={w}_{1}& 0& \cdots & 0\\ 0& {\gamma}_{2j}={w}_{2}& 0& 0\\ \vdots & \vdots & \ddots & \vdots \\ 0& 0& \cdots & {\gamma}_{mj}={w}_{m}\end{array}\right],

(19)

where, *w*_{
k
}∈ {0,1}, for *k* = 1, 2, ..., *m*. Because the prior probability of each *γ*_{
kj
}is independent, the joint prior probability for all possible formations can be written as p({\mathbf{\Phi}}_{j}={W}_{l})={\displaystyle \mathbf{\prod}_{k=1}^{m}p({\gamma}_{kj}={w}_{k})}. Then the conditional posterior probability of **Φ**_{
j
}can be written as

p({\mathbf{\Phi}}_{j}={W}_{j,{w}_{1}{w}_{2}\cdots {w}_{m}}|\cdots )=\frac{p({\mathbf{\Phi}}_{j}={W}_{j,{w}_{1}{w}_{2}\cdots {w}_{m}}){\displaystyle {\prod}_{i=1}^{n}f({y}_{i}|{\mathbf{\Phi}}_{j}={W}_{j,{w}_{1}{w}_{2}\cdots {w}_{m}},\cdots )}}{{\displaystyle {\sum}_{{g}_{1}\in \{0,1\}}{\displaystyle {\sum}_{{g}_{2}\in \{0,1\}}\cdots {\displaystyle {\sum}_{{g}_{m}\in \{0,1\}}\left(p({\mathbf{\Phi}}_{j}={W}_{j,{g}_{1}{g}_{2}\cdots {g}_{m}}){\displaystyle {\prod}_{i=1}^{n}f({y}_{i}|{\mathbf{\Phi}}_{j}={W}_{j,{g}_{1}{g}_{2}\cdots {g}_{m}},\cdots )}\right)}}}}.

(20)

The approach to sample **Φ**_{
j
}is similar to QTL genotypes sampling previously mentioned.

### Post-MCMC analysis

For summarizing the posterior sample, we use the mean of the posterior sample to estimate the QTL effect and the residual (co)variance, and the mode of the posterior probability or the peak of the 2log_{e}BF statistic to localize QTL. 2log_{e}BF statistic was introduced by Yi et al.[17] into QTL mapping, and BF statistic is defined as the ratio of the posterior odds to the prior odds for inclusion against exclusion of the locus [24]. The critical value of BF is 3 or 2log_{e}BF = 2.1 for declaring the existence of a QTL.

In single-trait analysis, we can pick the QTL by plotting the profile of the posterior probability or 2log_{e}BF statistic against the genome. In multitrait analysis, if only two traits are considered jointly, we can use a three-dimension graph to summarize the statistic for all traits jointly (e.g., Figure 2 in [19]). However, if the number of trait is greater than 2, we can't plot them in one graph. Instead, we can solve the problem by plotting the marginal posterior probability distribution. If we divide the genome into *H* bins, and denote each bin of *k* th trait with *ζ*_{
kg
}, for *g* = 1,2, ..., *H*, then the marginal posterior probability distribution of *ζ*_{
kg
}is defined as *p*(*ζ*_{
kg
}|**y**) = *p*[(*ζ*_{
kg
}= *λ*_{
kq
}) ∩ (*γ*_{
kq
}= 1)], where, *q* indicates the *q* th interval that locus *ζ*_{
kg
}resides in. Then \text{BF}({\zeta}_{kg})=\frac{p({\zeta}_{kg}|y)}{1-p({\zeta}_{kg}|y)}\cdot \frac{1-p({\zeta}_{kg})}{p({\zeta}_{kg})}, which can be calculated at each possible locus for each trait, respectively.