next up previous


This LaTeX document is available as postscript or asAdobe PDF.

Multiple Trait Animal Models
L. R. Schaeffer, March 1999

Animals are commonly observed for more than one trait because many traits affect overall profitability. Dairy cattle, for example, are observed for production traits (milk, fat, and protein yields), conformation traits (too many to list), calving ease, milking speed, temperament, survival, and disease susceptibility. Beef, swine, and sheep are observed for a number of weight traits, reproductive performance, litter size, carcass traits, and others. A multiple trait (MT) model is one in which two or more traits are analyzed simultaneously in order to take advantage of genetic and environmental correlations between traits.

Multiple trait models are useful for traits where the difference between genetic and residual correlations are large ( e.g. greater than .5 difference ) or where one trait has a much higher heritability than the other trait. In the latter case, traits with low heritability tend to gain more in accuracy than high heritability traits, although all traits benefit to some degree from the simultaneous analysis. Another use of MT models is for traits that occur at different times in the life of the animal, such that culling of animals results in fewer observations on animals for traits that occur later in life compared to those at the start. Consequently, animals which have observations later in life tend to have been selected based on their performance for earlier traits. Thus, analysis of later life traits by themselves could suffer from the effects of culling bias, and the resulting EBV could lead to errors in selecting future parents. An MT analysis that includes all observations on an animal upon which culling decisions have been based, has been shown to account for the selection that has taken place, and therefore gives unbiased estimates of breeding values for all traits.

MT models do not offer great increases in accuracy for cases where heritabilities of traits are similar in magnitude, and where both genetic and residual correlations are relatively the same. However, if culling bias exists, then an MT analysis should be performed even if the parameters are similar. If all animals are observed for all traits, then there would be no need to worry about culling bias because all traits would be equally affected.

An MT analysis relies on the accuracy of the genetic and residual correlations that are assumed. If the parameter estimates are greatly different from the underlying, unknown true values, then an MT analysis could do as much harm as it might do good.

Lastly, the researcher needs to consider the increased costs of computing MT analyses. Programs are more complicated, more memory and disk storage are usually needed, and verification of results might be more complicated. These have to be balanced against the benefits of an MT analysis. If culling bias is the main concern, then an MT model must be used regardless of the costs or no analysis should be done at all, except for the traits not affected by culling bias.

Models

Consider two traits with a single observation per trait on animals. A model should be specified separately for each trait. Usually, the same model is assumed for each trait, and this can greatly simplify the computational aspects, but such an assumption may be unrealistic in many situations.

Let the model equation for trait 1 be

y1ij = B1i + a1j + e1ij,

where B1i is a fixed effect with pB levels, a1j is a random, animal additive genetic effect for trait 1, and e1ij is a random residual environmental effect for trait 1. The model equation for trait 2 might be

y2ij = C2i + a2j + e2ij,

where C2i is a fixed effect (different from B1i for trait 1) with pC levels, a2j is a random, animal additive genetic effect for trait 2, and e2ij is a random residual environmental effect for trait 2.

For example, y1ij could be a trait like birthweight, so that B1i could identify animals born in the same season. Trait 2 could be yearling weights and C2i could identify contemporary groups of animals of the same sex, same herd, and same rearing unit within herd.

Because the two traits will be analyzed simultaneously, the variances and covariances need to be specified for the traits together. For example, the additive genetic variance-covariance (VCV) matrix could be written as

\begin{displaymath}{\bf G} = \left( \begin{array}{rr}
g_{11} & g_{12} \\ g_{12} ...
...
\left( \begin{array}{rr} 1 & 2 \\ 2 & 15 \end{array} \right), \end{displaymath}

and the residual environmental VCV matrix as

\begin{displaymath}{\bf R} = \left( \begin{array}{rr}
r_{11} & r_{12} \\ r_{12}...
...left( \begin{array}{rr} 10 & 5 \\ 5 & 100 \end{array} \right). \end{displaymath}

The genetic and residual correlations are, respectively,

\begin{eqnarray*}\rho_{g} & = & 2 / (15)^{.5} = .516, \\
\rho_{r} & = & 5 / (1000)^{.5} = .158
\end{eqnarray*}


with

\begin{displaymath}h^{2}_{1} = \frac{1}{11} = .0909, \end{displaymath}

and

\begin{displaymath}h^{2}_{2} = \frac{15}{115} = .1304. \end{displaymath}

For all data, then

\begin{displaymath}Var \left( \begin{array}{c} {\bf a}_{1} \\ {\bf a}_{2}
\end{...
...A}g_{12} \\ {\bf A}g_{12} & {\bf A}g_{22}
\end{array} \right). \end{displaymath}

The structure of the residual VCV matrix over all observations can be written several ways depending on whether allowance is made for missing observations on either trait for some animals. If all animals were observed for both traits, then

\begin{displaymath}Var \left( \begin{array}{c} {\bf e}_{1} \\ {\bf e}_{2}
\end{...
...}r_{12} \\ {\bf I}r_{12} & {\bf I}r_{22}
\end{array} \right). \end{displaymath}

Simulating Records

When simulating data for a multiple trait problem it is best to generate observations for all animals for all traits. Then one can go through the simulated data and randomly delete observations to simulate a missing data situation, or selectively delete observations to imitate culling decisions. Another simplification is to assume that the model for each trait is the same, and then for a factor that does not belong with a given trait just make the true values of levels for that factor and trait equal to zero. In matrix form, the model equation for one animal would be

\begin{eqnarray*}\left( \begin{array}{c} y_{1ij} \\ y_{2ij} \end{array} \right) ...
...{array}{c} \mbox{ RND(3)} \\ \mbox{ RND(4)}
\end{array} \right),
\end{eqnarray*}


where ${\bf L}_{G}$ and ${\bf L}_{R}$ are lower triangular matrices such that

\begin{eqnarray*}{\bf G} & = & {\bf L}_{G}{\bf L}'_{G} \\
& = & \left( \begin{...
... 0 \\ (2.5)^{.5} &
(97.5)^{.5} \end{array} \right){\bf L}'_{R}.
\end{eqnarray*}


Let B11 = 6.7 and B12 = 6.3 for trait 1, and because factor B is not in the model for trait 2, then B21=0 and B22=0. Similarly, C21=25, C22=40, and C23=55 for trait 2, and because factor C is not in the model for trait 1, then C11=0, C12=0, and C13=0. Suppose the animal is a base animal, then bii=1 and the parent averages for traits 1 and 2 are assumed to be zero, then the observations would be

\begin{eqnarray*}\left( \begin{array}{c} y_{1ij} \\ y_{2ij} \end{array} \right) ...
...& \left( \begin{array}{r} 5.7106 \\ 26.8771 \end{array} \right).
\end{eqnarray*}


The following data (rounded off) were simulated according to the preceeding scheme and parameters.

Animal Sire Dam B-level C-level Trait 1 Trait 2
1 0 0 1 1 2.3 39
2 0 0 1 2 2.6 39
3 0 0 1 3 9.8 53
4 0 0 1 1 4.7 4
5 0 0 1 2 5.5 63
6 1 3 2 3 2.5 64
7 1 4 2 2 8.4 35
8 1 5 2 3 8.2 41
9 2 3 2 1 9.0 27
10 2 4 2 1 7.8 32
11 2 5 2 2 2.8 46
12 6 10 2 3 7.4 67

To simulate selection, assume that all animals had trait 1 observed, but for any animal with a trait 1 value below 3.0, then their trait 2 observation was removed. Four trait 2 observations were deleted, giving the results in the table below.

Animal Sire Dam B-level C-level Trait 1 Trait 2
1 0 0 1 1 2.3  
2 0 0 1 2 2.6  
3 0 0 1 3 9.8 53
4 0 0 1 1 4.7 4
5 0 0 1 2 5.5 63
6 1 3 2 3 2.5  
7 1 4 2 2 8.4 35
8 1 5 2 3 8.2 41
9 2 3 2 1 9.0 27
10 2 4 2 1 7.8 32
11 2 5 2 2 2.8  
12 6 10 2 3 7.4 67

HMME

Organize the data by traits within animals. With two traits there are three possible residual matrices per animal, i.e.,

\begin{eqnarray*}{\bf R}_{12} & = & \left( \begin{array}{rr} 10 & 5 \\ 5 & 100
...
... \left( \begin{array}{rr} 0 & 0 \\ 0 & 100
\end{array} \right),
\end{eqnarray*}


depending on whether both traits, trait 1, or trait 2, respectively, were observed. In the example data, only ${\bf R}_{12}$ and ${\bf R}_{1}$ are needed. To simplify notation, let

\begin{eqnarray*}{\bf E}_{12} & = & {\bf R}^{-1}_{12} \\
& = & \frac{1}{975} \...
... .102564 & -.005128 \\
-.005128 & .010256 \end{array} \right),
\end{eqnarray*}


and

\begin{eqnarray*}{\bf E}_{1} & = & {\bf R}^{-}_{1} \\
& = & \left( \begin{arra...
... \left( \begin{array}{rr} 0 & 0 \\ 0 & .01
\end{array} \right).
\end{eqnarray*}


Again, to simplify construction of the MME, pretend that both traits have the same model equation, so that

ytijk = Bti+Ctj+atk+etijk.

There are 2 levels of factor B, three levels of factor C, and 12 animals. For a single trait model this would give MME of order 17. Construct a table of order 17. The elements of this table will be matrices of order 2(in general, order t). Start with animal 1, then

\begin{displaymath}{\bf y}_{1} = \left( \begin{array}{c} 2.3 \\ - \end{array}
...
..._{1} =
\left( \begin{array}{c} .23 \\ 0 \end{array} \right). \end{displaymath}

Now add ${\bf E}_{1}$ to the boxes in the MME table as follows:

  B1 B2 C1 C2 C3 a1 $\cdots$ RHS
B1 ${\bf E}_{1}$   ${\bf E}_{1}$     ${\bf E}_{1}$   ${\bf E}_{1}{\bf y}_{1}$
B2                
C1 ${\bf E}_{1}$   ${\bf E}_{1}$     ${\bf E}_{1}$   ${\bf E}_{1}{\bf y}_{1}$
C2                
C3                
a1 ${\bf E}_{1}$   ${\bf E}_{1}$     ${\bf E}_{1}$   ${\bf E}_{1}{\bf y}_{1}$

Similarly for animal 2,

\begin{displaymath}{\bf y}_{2} = \left( \begin{array}{c} 2.6 \\ - \end{array}
...
..._{2} =
\left( \begin{array}{c} .26 \\ 0 \end{array} \right). \end{displaymath}

Accumulating into the MME table gives

  B1 B2 C1 C2 C3 $\cdots$ a2 RHS
B1 $2{\bf E}_{1}$   ${\bf E}_{1}$ ${\bf E}_{1}$     ${\bf E}_{1}$ ${\bf E}_{1}({\bf y}_{1}+{\bf y}_{2})$
B2                
C1 ${\bf E}_{1}$   ${\bf E}_{1}$         ${\bf E}_{1}{\bf y}_{1}$
C2 ${\bf E}_{1}$     ${\bf E}_{1}$     ${\bf E}_{1}$ ${\bf E}_{1}{\bf y}_{2}$
C3                
a2 ${\bf E}_{1}$     ${\bf E}_{1}$     ${\bf E}_{1}$ ${\bf E}_{1}{\bf y}_{2}$

For animal 3,

\begin{displaymath}{\bf y}_{3} = \left( \begin{array}{c} 9.8 \\ 53 \end{array}
...
...
\left( \begin{array}{c} .7333 \\ .4933 \end{array} \right), \end{displaymath}

and the MME table becomes

  B1 B2 C1 C2 C3 $\cdots$ a3 RHS
B1 $2{\bf E}_{1}+{\bf E}_{12}$   ${\bf E}_{1}$ ${\bf E}_{1}$ ${\bf E}_{12}$   ${\bf E}_{12}$ ${\bf E}_{1}({\bf y}_{1}+{\bf y}_{2})+{\bf E}_{12}{\bf y}_{3}$
B2                
C1 ${\bf E}_{1}$   ${\bf E}_{1}$         ${\bf E}_{1}{\bf y}_{1}$
C2 ${\bf E}_{1}$     ${\bf E}_{1}$       ${\bf E}_{1}{\bf y}_{2}$
C3 ${\bf E}_{12}$       ${\bf E}_{12}$   ${\bf E}_{12}$ ${\bf E}_{12}{\bf y}_{3}$
a3 ${\bf E}_{12}$       ${\bf E}_{12}$   ${\bf E}_{12}$ ${\bf E}_{12}{\bf y}_{3}$

The remaining animals are processed in the same manner. The resulting equations are of order 34 by 34. To these ${\bf A}^{-1}\otimes
{\bf G}^{-1}$ must be added to the animal by animal submatrix in order to form the full HMME. However, solutions for the B-factor for trait 2 are not needed because the B-factor does not affect trait 2, and solutions for the C-factor for trait 1 are not needed because the C-factor does not affect trait 1. Therefore, remove rows (and columns) 2, 4, 5, 7, and 9, or if an iterative solution is being computed, then require that the solutions for B21, B22, C11, C12, and C13 are always equal to zero. The solutions to the HMME, for this example, were

\begin{eqnarray*}B_{11} & = & 5.0209 \\
B_{12} & = & 6.5592 \\
C_{21} & = & 20.0882 \\
C_{22} & = & 49.0575 \\
C_{23} & = & 51.9553
\end{eqnarray*}


The animal additive genetic solutions are shown in the table below.

Animal Sire Dam Trait 1 Trait 2
1 0 0 -.3573 -1.6772
2 0 0 -.0730 1.0418
3 0 0 .4105 1.1707
4 0 0 -.0449 -1.4922
5 0 0 .0646 .9570
6 1 3 -.1033 -.1410
7 1 4 -.1975 -2.2983
8 1 5 -.1410 -.9633
9 2 3 .3079 1.6227
10 2 4 .1426 1.1273
11 2 5 -.1830 .6418
12 6 10 .1554 1.5089

The correlation between the animal additive genetic solutions for traits 1 and 2 was .74 which is greater than the .52 assumed in the original ${\bf G}$.

Partitioning the MT solution

Some insight into the workings of multiple trait HMME can be gained by partitioning an animal's additive genetic solution for any one trait. Take animal 10 from the previous example, because this animal has a record on both traits and progeny, plus parents are known. The partitioning results in contributions attributable to each trait weighted by the genetic correlations between traits. In addition to the Data, Parent Average, and Progeny contributions, there is also a contribution from the direct genetic solutions for the other traits. The partitions will be presented in tabular form.

Partitions of Solution for Animal 10
Contribution Trait 1 Trait 2
Data $r^{11}(y_{10,1}-\hat{B}_{12})$ $r^{12}(y_{10,2}-\hat{C}_{21})$
  (.102564)(7.8-6.5592) (-.005128)(32-20.0882)
  .12726141 -.06108371
Parent $2g^{11}(.5)(\hat{a}_{2,1}+\hat{a}_{4,1})$ $2g^{12}(.5)(\hat{a}_{2,2}+\hat{a}_{4,2})$
Average 2.72727(-.0589225) -.363636(-.2252442)
  -.16069773 .08190698
Progeny $.5g^{11}(2\hat{a}_{12,1}-\hat{a}_{6,1})$ $.5g^{12}(2\hat{a}_{12,2}-\hat{a}_{6,2})$
  .681818(.41411480) -.090909(3.15877680)
  .282351 -.28716153
Other   $-(2.5g^{12}+r^{12})\hat{a}_{10,2}$
Traits   - (-.459674)(1.1272952)
    .51818829

Summing all of the above pieces gives $(2.5g^{11}+r^{11})\hat{a}_{10,1}$. There would be additional columns if more traits were included in the analysis, i.e. one for each additional trait. One could interpret the results in the above table as follows. The Data contribution from trait 1 to the trait 1 solution was .1273, but the trait 2 correlated information said that it should be .0611 lower. The parent average for trait 1 contributed -.1607, but the correlated parent average from trait 2 indicated it should be .0819 higher. The Progeny contribution from trait 1 was nearly equal and opposite to the contribution from trait 2. The contribution from the direct genetic solution for trait 2 was seemingly large.

Another way to look at the table is to combine all of the partitions of the trait 2 column into one figure, because the Data, Parent Average, and Progeny contributions from trait 2 are components of the additive genetic solution for trait 2, but they are weighted slightly differently, that is,

\begin{displaymath}(2.5 g^{12}+r^{12})\hat{a}_{10,2} \approx
r^{12}(y_{10,2}-\...
...2}
+ \hat{a}_{4,2})+ .5g^{12}(2\hat{a}_{12,2}-\hat{a}_{6,2}). \end{displaymath}

The difference between the two (which is the total of that column) should be fairly small in most situations, but does contribute towards the trait 1 additive genetic solution in a regressed fashion. The weights on the partitions in the trait 1 column sum to 1, but the weights on the partitions in the trait 2 column sum to 0. Bruce Tier's MT Tricks

There is another way to construct the MME without the need of forming different inverses of ${\bf R}$ for missing traits. If a trait is missing, then that observation is assigned to its own contemporary group in the model for that trait. In the example data there were four missing observations. Animal 1 would be assigned to C24, animal 2 to C25, animal 6 to C26 and animal 11 to C27, respectively. In this case only trait 2 observations were missing. If trait 1 observations were also missing, then animals would be assigned to separate levels of factor B. In this way, only one residual VCV matrix is needed, i.e. ${\bf R}_{12}$. Let ${\bf X}$ represent the design matrix for fixed effects (factors B and C) for either trait. Note the four extra columns for factor C for the animals with missing trait 2 observations.

\begin{displaymath}{\bf Xb} = \left( \begin{array}{rrrrrrrrr}
1 & 0 & 0 & 0 & 0 ...
... \\
C_{24} \\ C_{25} \\ C_{26} \\ C_{27} \end{array} \right). \end{displaymath}

A missing observation is replaced with zero. The resulting solutions are identical to those given earlier, except that now there are also solutions for the single observation contemporary groups. That is, C24=2.8590, C25=.1322, C26=2.1189, and C27=1.1463. These solutions do not mean anything, but must be calculated.

To prove that this trick will work, take ${\bf R}^{-1}$ and do a Gaussian elimination (i.e. absorption) of the row and column corresponding to the missing trait, say trait 2,

\begin{displaymath}\left( \begin{array}{rr} e^{11} & e^{12} \\ e^{12} & e^{22}
...
... \left( \begin{array}{rr}
e^{12} & e^{22} \end{array} \right), \end{displaymath}


\begin{displaymath}= \left( \begin{array}{rr} e^{11} & e^{12} \\ e^{12} & e^{22}...
...^{-1}e^{12} & e^{12} \\
e^{12} & e^{22} \end{array} \right), \end{displaymath}


\begin{displaymath}= \left( \begin{array}{cc} e^{11}-e^{12}(e^{22})^{-1}e^{12} & 0 \\
0 & 0 \end{array} \right). \end{displaymath}

Recall that for a matrix of order 2 that

\begin{eqnarray*}e^{11} & = & e_{22}/\mid {\bf R} \mid, \\
e^{12} & = & -e_{12}...
...R} \mid, \\
\mid {\bf R} \mid & = & (e_{11}e_{22}-e_{12}e_{12})
\end{eqnarray*}


then

\begin{eqnarray*}e^{11}-e^{12}(e^{22})^{-1}e^{12} & = &
( e_{22} - e_{12}(e_{1...
...12}/e_{11}(e_{11}e_{22}-
e_{12}e_{12}) \\
& = & (e_{11})^{-1}
\end{eqnarray*}


which is exactly the weight applied to records on animals with only trait 1 observed. This proof can be extended to any number of traits recorded and any number missing, by partitioning ${\bf R}$ into

\begin{displaymath}\left( \begin{array}{cc} {\bf R}_{oo} & {\bf R}_{om} \\
{\bf R}_{mo} & {\bf R}_{mm} \end{array} \right), \end{displaymath}

where the subscript o refers to traits that were observed and m refers to traits that were missing on an animal. Then it can be easily shown that

\begin{displaymath}{\bf R}^{-1}_{oo} = {\bf R}^{oo}-{\bf R}^{om}({\bf R}^{mm})^{-1}
{\bf R}^{mo}. \end{displaymath}

This trick is not very practical, for example, when one trait has 1 million observations and trait 2 has only 100,000 observations, then there would be 900,000 extra single observation subclasses created for trait 2. However, if the percentages of missing observations are relatively small, or if many traits are being considered, then pretending all observations are present may make programming easier.

Estimation of Covariances

Derivative free REML is one option for estimating variances and covariances in a multi-trait situation. The EM algorithm is not suitable due to the requirement for the traces of inverse elements that are needed. Even DF REML takes considerably more time as the number of parameters to be estimated increases.

Another option is the Bayesian approach, where operations are performed in t dimensions, for t being the number of traits. Thus, for a solution to the MME, the $t \times 1$ vector for any one fixed effect, for example, would be

\begin{displaymath}\hat{\beta}_{i} = ({\bf X}'_{i}{\bf R}^{-1}{\bf X}_{i})^{-1}
...
...f X}'_{i}{\bf R}^{-1}( {\bf y}_{i} - {\bf W}_{-i}\beta_{-i} ), \end{displaymath}

then

\begin{displaymath}\beta_{i} = \hat{\beta}_{i} + {\bf L}{\bf v}, \end{displaymath}

where

\begin{displaymath}{\bf LL}' = ({\bf X}'_{i}{\bf R}^{-1}{\bf X}_{i})^{-1}, \end{displaymath}

and ${\bf v}$ is a $t \times 1$ vector of random normal deviates. Similar formulas can be derived for the random factors of the model. The conditional distributions for these factors are assumed to be normally distributed.

If ${\bf a}_{i}$ is the $q \times 1$ vector of animal solutions for trait i, then form

\begin{displaymath}{\bf U} = \left( \begin{array}{rrrr}
{\bf a}_{1} & {\bf a}_{2} & \cdots & {\bf a}_{t} \end{array} \right), \end{displaymath}

followed by

\begin{displaymath}{\bf S}_{a} = ({\bf U}'{\bf A}^{-1}{\bf U} + \nu_{a} {\bf G}_{a}), \end{displaymath}

which is then inverted and a Cholesky decomposition is applied to the inverse, i.e.

\begin{displaymath}{\bf L}_{a}{\bf L}'_{a} = {\bf S}^{-1}_{a}, \end{displaymath}

where ${\bf L}_{a}$ is supplied to a Wishart distribution random number generator to give a new sample matrix for the inverse of the additive genetic variances and covariances.

A difficult part of a multiple trait analysis, when missing traits are possible, is the calculation of the appropriate residual matrix of sums of squares and cross products. The residual effect for any one trait is

\begin{displaymath}e_{ti} = y_{ti} - {\bf w}'_{i}\beta_{t}, \end{displaymath}

and for a particular animal, k,

\begin{displaymath}{\bf res}_{k} = \left( \begin{array}{cc} {\bf e}'_{o} & {\bf e}'_{m}
\end{array} \right), \end{displaymath}

where the subscripts o and m refer to observed and missing, respectively, and where ${\bf e}'_{m}={\bf0}'$. In order to calculate the correct sum of squares and crossproducts of the residuals for REML or the Bayesian approach, a prediction of ${\bf e}_{m}$ is needed. This can be easily done by

\begin{eqnarray*}{\bf res}'_{k*} & = & \left( \begin{array}{ll} {\bf R}_{oo} &
...
...ft( \begin{array}{c}
{\bf e}_{o} \\ {\bf0} \end{array} \right).
\end{eqnarray*}


For the first animal in the example data, for example, the estimated residual for trait 1 was -2.363638, then

\begin{eqnarray*}{\bf res}'_{1*} & = & \left( \begin{array}{rr} 10 & 5 \\ 5 & 10...
...t( \begin{array}{r} -2.363638 \\ -1.181819
\end{array} \right).
\end{eqnarray*}


If you use Bruce Tier's MME with the single observation contemporary groups for missing trait observations, then the residuals can be calculated directly by $({\bf y}-{\bf W}\beta)$ using zero as the observation for the missing traits and using the solutions for the single observation contemporary groups. This gives the exact same residual estimates as the above methodology. Therefore, Tier's approach is handy for the Gibb's sampling algorithm.

Once the residuals are calculated for all animals with records, then calculate

\begin{displaymath}{\bf res}'_{k*}{\bf res}_{k*}, \end{displaymath}

and sum across N animals. Then let

\begin{displaymath}{\bf S}_{e} = ( \sum_{k=1}^{N}({\bf res}'_{k*}{\bf res}_{k*})
+ \nu_{e} {\bf R}_{e}), \end{displaymath}

which is then inverted and a Cholesky decomposition is applied to the inverse, i.e.

\begin{displaymath}{\bf L}_{e}{\bf L}'_{e} = {\bf S}^{-1}_{e}, \end{displaymath}

where ${\bf L}_{e}$ is supplied to a Wishart distribution random number generator to give a new sample matrix for the inverse of the residual variances and covariances.

Simplifications

The previous description applies to general multiple trait models where the model for each trait can be different and where animals need not be observed for each trait. This is frequently the true situation when a multiple trait analysis is utilized in order to account for selection bias. However, there are special situations where multiple trait analyses can be greatly simplified in terms of the necessary calculations. One special case is when all traits have the exact same model and when all animals are observed for all traits, and only additive genetic and residual effects are random in the model. In this case a canonical transformation can be applied to the traits. That is, a matrix ${\bf M}$ can be found such that ${\bf M}'{\bf R}{\bf M}={\bf I}$ and that ${\bf M}'{\bf G}{\bf M}$ is a diagonal matrix. Thus, the transformed traits are considered to be uncorrelated, and consequently each transformed trait can be analyzed independently of the others. At the end, the results can be back-transformed to the original scale using ${\bf M}^{-1}$. The steps to find ${\bf M}$ are as follows.

Step 1.
Find the eigenvalues and eigenvectors of ${\bf R}$.

\begin{displaymath}{\bf R} = {\bf UBU}', \end{displaymath}

such that ${\bf U}'{\bf U} = {\bf I}$ and ${\bf B}$ is diagonal with all diagonals greater than 0.
Step 2.
Form ${\bf P}={\bf B}^{-.5}{\bf U}'$, then

\begin{eqnarray*}{\bf PRP}' & = & {\bf B}^{-.5}{\bf U}'({\bf UBU}'){\bf UB}^{-.5...
...-.5}{\bf B}^{.5})({\bf B}^{-.5}{\bf B}^{.5}) \\
& = & {\bf I}
\end{eqnarray*}


Step 3.
Apply ${\bf P}$ to ${\bf G}$, and find the eigenvalues and eigenvectors of that product.

\begin{displaymath}{\bf PGP}' = {\bf LDL}', \end{displaymath}

such that ${\bf L}'{\bf L}={\bf I}$ and ${\bf D}$ is diagonal with all diagonal elements greater than zero.
Step 4.
Form ${\bf M}={\bf L}'{\bf P}$, then

\begin{eqnarray*}{\bf MRM}' & = & {\bf L}'{\bf PRP}'{\bf L} \\
& = & {\bf L}'{...
...& ({\bf L}'{\bf L}){\bf D}({\bf L}'{\bf L}) \\
& = & {\bf D}.
\end{eqnarray*}


The data on animal k are transformed as ${\bf My}_{k} = {\bf z}_{k}$, and each transformed trait is analyzed separately. An estimate of the residual variance for any transformed trait should be close to 1, and the estimate of the additive genetic variance should be close to the corresponding diagonal element of ${\bf D}$. Transform the estimates back to the original scale by premultiplying by ${\bf M}^{-1}$ and postmultiplying by ${\bf M}^{-T}$. This process needs to be continued iteratively until convergence is reached.

The canonical transformation cannot be applied when the model contains more than additive genetic and residual variance-covariance matrices because it is impossible to diagonalize more than two matrices. However, there are approximation methods for diagonalizing more than two matrices at a time which have been applied to multiple trait situations. These methods are not covered in these notes.


next up previous

This LaTeX document is available as postscript or asAdobe PDF.

Larry Schaeffer
1999-03-29