This LaTeX document is available as postscript or asAdobe PDF.

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

where

where

For example, *y*_{1ij} could be a trait like birthweight, so that
*B*_{1i} could identify animals born in the same season. Trait 2
could be yearling weights and *C*_{2i} 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

and the residual environmental VCV matrix as

The genetic and residual correlations are, respectively,

with

and

For all data, then

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

**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

where and are lower triangular matrices such that

Let
*B*_{11} = 6.7 and
*B*_{12} = 6.3 for trait 1, and because factor
B is not in the model for trait 2, then *B*_{21}=0 and *B*_{22}=0.
Similarly, *C*_{21}=25, *C*_{22}=40, and *C*_{23}=55 for trait 2, and
because factor C is not in the model for trait 1, then
*C*_{11}=0, *C*_{12}=0, and *C*_{13}=0. Suppose the animal is a
base animal, then *b*_{ii}=1 and the parent averages for traits 1 and
2 are assumed to be zero, then the observations would be

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.,

depending on whether both traits, trait 1, or trait 2, respectively, were observed. In the example data, only and are needed. To simplify notation, let

and

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

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

Now add to the boxes in the MME table as follows:

B_{1} |
B_{2} |
C_{1} |
C_{2} |
C_{3} |
a_{1} |
RHS | ||

B_{1} |
||||||||

B_{2} |
||||||||

C_{1} |
||||||||

C_{2} |
||||||||

C_{3} |
||||||||

a_{1} |

Similarly for animal 2,

Accumulating into the MME table gives

B_{1} |
B_{2} |
C_{1} |
C_{2} |
C_{3} |
a_{2} |
RHS | ||

B_{1} |
||||||||

B_{2} |
||||||||

C_{1} |
||||||||

C_{2} |
||||||||

C_{3} |
||||||||

a_{2} |

For animal 3,

and the MME table becomes

B_{1} |
B_{2} |
C_{1} |
C_{2} |
C_{3} |
a_{3} |
RHS | ||

B_{1} |
||||||||

B_{2} |
||||||||

C_{1} |
||||||||

C_{2} |
||||||||

C_{3} |
||||||||

a_{3} |

The remaining animals are processed in the same manner. The resulting
equations are of order 34 by 34. To these
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 *B*_{21},
*B*_{22}, *C*_{11}, *C*_{12}, and *C*_{13} are always equal to zero.
The solutions to the HMME, for this example, were

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 .

**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 | ||

(.102564)(7.8-6.5592) | (-.005128)(32-20.0882) | |

.12726141 | -.06108371 | |

Parent | ||

Average | 2.72727(-.0589225) | -.363636(-.2252442) |

-.16069773 | .08190698 | |

Progeny | ||

.681818(.41411480) | -.090909(3.15877680) | |

.282351 | -.28716153 | |

Other | ||

Traits | - (-.459674)(1.1272952) | |

.51818829 |

Summing all of the above pieces gives . 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,

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.

There is another way to construct the MME without the need of
forming different inverses of
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 *C*_{24},
animal 2 to *C*_{25}, animal 6 to *C*_{26} and animal 11 to *C*_{27},
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.
.
Let
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.

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,

To prove that this trick will work, take
and do
a Gaussian elimination (i.e. absorption) of the row and column
corresponding to the missing trait, say trait 2,

Recall that for a matrix of order 2 that

then

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 into

where the subscript

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
vector for any one fixed effect,
for example, would be

then

where

and is a 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
is the
vector of animal solutions for
trait *i*, then form

followed by

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

where 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

and for a particular animal,

where the subscripts

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

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 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

and sum across

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

where 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
can be
found such that
and that
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
.
The steps to find
are as
follows.

**Step 1.**- Find the eigenvalues and eigenvectors of .

such that and is diagonal with all diagonals greater than 0. **Step 2.**- Form
,
then

**Step 3.**- Apply
to ,
and find the eigenvalues
and eigenvectors of that product.

such that and is diagonal with all diagonal elements greater than zero. **Step 4.**- Form
,
then

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.

This LaTeX document is available as postscript or asAdobe PDF.