next up previous


This LaTeX document is available as postscript or asAdobe PDF.

Maternal Effects Model
L. R. Schaeffer, March 1999
Updated March 16, 2000

In some species of livestock, such as beef cattle, sheep or swine, the female provides an environment for its offspring to survive and grow. Females vary in their ability to provide a good environment for their offspring, and this variability has a genetic basis. The offspring inherit directly an ability to grow (or survive) from both parents, and environmentally do better or poorer depending on their dam's maternal ability. Maternal ability is a genetic trait and is transmitted, as usual, from both parents, but maternal ability is only expressed by females when they have a calf (i.e. much like milk yield in dairy cows).

A model to account for maternal ability is

\begin{displaymath}{\bf y} = {\bf Xb} + {\bf Z}_{1}{\bf a} + {\bf Z}_{2}{\bf m} +
{\bf Z}_{3}{\bf p} + {\bf e}, \end{displaymath}

where ${\bf y}$ is the growth trait of a young animal, ${\bf b}$ is a vector of fixed factors influencing growth, such as contemporary group, sex of the offspring, or age of dam, ${\bf a}$ is a vector of random additive genetic effects (i.e. direct genetic effects) of the animals, ${\bf m}$ is a vector of random maternal genetic (dam) effects, and ${\bf p}$, in this model, is a vector of maternal permanent environmental effects (because dams may have more than one offspring in the data).

The expectations of the random vectors, ${\bf a}$, ${\bf m}$, ${\bf p}$, and ${\bf e}$ are all null vectors in a model without selection, and the variance-covariance structure is

\begin{displaymath}Var \left( \begin{array}{c} {\bf a} \\ {\bf m} \\ {\bf p} \\ ...
... {\bf0} & {\bf0} & {\bf I}\sigma^{2}_{e}
\end{array} \right), \end{displaymath}

where $\sigma^{2}_{a}$ is the additive genetic variance, $\sigma^{2}_{m}$ is the maternal genetic variance, $\sigma_{am}$ is the additive genetic by maternal genetic covariance, and $\sigma^{2}_{p}$ is the maternal permanent environmental variance. Also,

\begin{displaymath}\left( \begin{array}{c\vert l} {\bf a} & \\
& {\bf A},{\bf ...
...rray} \right), &
{\bf G} \otimes {\bf A} \end{array} \right), \end{displaymath}

where

\begin{displaymath}{\bf G} = \left( \begin{array}{cc}
\sigma^{2}_{a} & \sigma_{am} \\ \sigma_{am} & \sigma^{2}_{m}
\end{array} \right), \end{displaymath}

and

\begin{displaymath}{\bf p} \mid {\bf I}, \sigma^{2}_{p} \sim N( {\bf0},
{\bf I}\sigma^{2}_{p} ), \end{displaymath}

and

\begin{displaymath}{\bf e} \sim N({\bf0},{\bf I}\sigma^{2}_{e}). \end{displaymath}

In this model, a female animal, i, could have its own growth record for estimating $\hat{a}_{i}$. The same female could later have offspring of its own for estimating $\hat{m}_{i}$ and $\hat{p}_{i}$, and the offspring would also contribute towards $\hat{a}_{i}$. The maternal effects model can be more complicated if, for example, embryo transfer is practiced. Recipient dams would have maternal effects, but would not have direct genetic effects on that calf, see Schaeffer and Kennedy (1989). To better understand this model we will simulate some records. Let

\begin{displaymath}{\bf G} = \left( \begin{array}{cc}
\sigma^{2}_{a} & \sigma_{...
...ft( \begin{array}{rr}
49 & -7 \\ -7 & 26 \end{array} \right). \end{displaymath}

Now any positive definite matrix can be partitioned into the product of a matrix times its transpose (i.e. Cholesky decomposition), or

\begin{eqnarray*}{\bf G} & = & {\bf L}{\bf L}' \\
{\bf L} & = & \left( \begin{array}{rr}
7 & 0 \\ -1 & 5 \end{array} \right).
\end{eqnarray*}


Let $\sigma^{2}_{p} = 9$ and $\sigma^{2}_{e} = 81 $. This model differs from previous models in that both the additive genetic and maternal genetic effects need to be generated simultaneously because these effects are genetically correlated. We will generate true breeding values for three animals, A, B, and C, where C is an offspring of sire A and dam B, and then we will generate an observation on animal C.

For A, we need two random normal deviates which will be pre-multiplied by ${\bf L}$. We assume that A is a base population animal that is unrelated to B. Let ${\bf w}' = ( 2.533 \ \ -.299 )$, then

\begin{eqnarray*}\left( \begin{array}{c} a_{A} \\ m_{A} \end{array} \right)
& =...
...& \left( \begin{array}{r} 17.731 \\ -4.028 \end{array} \right).
\end{eqnarray*}


Similarly for animal B,

\begin{eqnarray*}\left( \begin{array}{c} a_{B} \\ m_{B} \end{array} \right)
& ...
... & \left( \begin{array}{r} -7.987 \\ 2.316 \end{array} \right).
\end{eqnarray*}


Creating a progeny's true breeding value is similar to the scalar version. Take the average of the parents' true breeding values and add a random Mendelian sampling term.

\begin{eqnarray*}\left( \begin{array}{c} a_{C} \\ m_{C} \end{array} \right)
& =...
... & \left( \begin{array}{r}
6.233 \\ .371 \end{array} \right).
\end{eqnarray*}


Note that all animals have both a direct and maternal genetic breeding values.

Next for all dams a maternal permanent environmental effect should be generated. In this case only for animal B, multiply a random normal deviate by $\sigma_{p}= 3$, suppose it is -4.491.

An observation for animal C is created by following the model equation,

\begin{eqnarray*}y & = & \mbox{Fixed Effects} + a_{C} + m_{B} + p_{B} + \sigma_{...
... 140 + 6.233 + 2.316 + (-4.491) + (9)(1.074) \\
& = & 153.724.
\end{eqnarray*}


The Fixed Effects contribution of 140 was arbitrarily chosen for this example. The point of main importance is that the observation on animal C consists of the direct genetic effect of animal Cplus the maternal genetic effect of the dam (B) plus the maternal permanent environmental effect of the dam (B) plus a residual.

HMME

To illustrate the calculations, assume the data as given in the table below.

Animal Sire Dam CG Weight
5 1 3 1 156
6 2 3 1 124
7 1 4 1 135
8 2 4 2 163
9 1 3 2 149
10 2 4 2 138

CG stands for contemporary group, the only fixed effect in this example. Assume that the appropriate variance parameters are those which were used in the simulation in the previous section. Based on the matrix formulation of the model, the MME are

\begin{displaymath}\left( \begin{array}{cccc} {\bf X}'{\bf X} & {\bf X}'{\bf Z}_...
...bf Z}'_{2}{\bf y} \\ {\bf Z}'_{3}{\bf y}
\end{array} \right), \end{displaymath}

where

\begin{eqnarray*}\left( \begin{array}{cc} k_{11} & k_{12} \\ k_{12} & k_{22}
\e...
...rray}{rr} 1.7192 & .4628 \\ .4628 & 3.2400
\end{array} \right).
\end{eqnarray*}


Note that these numbers are not equal to

\begin{displaymath}\left( \begin{array}{rr} 81/49 & 81/(-7) \\ 81/(-7) & 81/26
\end{array} \right). \end{displaymath}

Finally, $k_{33} = \sigma^{2}_{e} / \sigma^{2}_{p} = 81/9 = 9$.

The matrices are

\begin{displaymath}{\bf X} = \left( \begin{array}{rr}
1 & 0 \\ 1 & 0 \\ 1 & 0 \...
...bf y} = \left( \begin{array}{r} 415 \\ 450 \end{array}\right), \end{displaymath}


\begin{displaymath}{\bf Z}_{1} = \left( \begin{array}{rrrrrrrrrr}
0 & 0 & 0 & 0...
...
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1
\end{array} \right), \end{displaymath}


\begin{displaymath}{\bf Z}_{2} = \left( \begin{array}{rrrrrrrrrr}
0 & 0 & 1 & 0...
...
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0
\end{array} \right), \end{displaymath}


\begin{displaymath}{\bf Z}_{3} = \left( \begin{array}{rr}
1 & 0 \\ 1 & 0 \\ 0 &...
...f y} = \left( \begin{array}{r} 429 \\ 436 \end{array} \right). \end{displaymath}

The other two right hand side matrices can be easily obtained from ${\bf y}$ and ${\bf Z}'_{3}{\bf y}$. Thus, the order of the MME will be 24. The inverse of the relationship matrix is

\begin{displaymath}{\bf A}^{-1} = \frac{1}{2} \left( \begin{array}{rrrrrrrrrr}
...
...
0 & -2 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 4 \end{array} \right). \end{displaymath}

The solutions to the MME are

\begin{displaymath}\hat{\bf b} = \left( \begin{array}{r} 137.8469 \\ 150.4864
\...
...\left( \begin{array}{r} .0658 \\ -.0658
\end{array} \right), \end{displaymath}


\begin{displaymath}\hat{\bf a} = \left( \begin{array}{r} 2.3295 \\ -2.3295 \\
...
...92 \\ -.1254 \\ -.3795 \\ .0136 \\ .4499
\end{array} \right). \end{displaymath}

No correlations with true values were calculated for this small example.

Estimation of Variances

REML

Let the inverse of HMME be represented as

\begin{displaymath}\left( \begin{array}{cccc}
{\bf C}_{bb} & {\bf C}_{ba} & {\b...
...bf C}_{pa} & {\bf C}_{pm} & {\bf C}_{pp}
\end{array} \right). \end{displaymath}

Then the REML formulas and results for this model are

\begin{eqnarray*}\hat{\sigma}^{2}_{e} & = & {\bf y}'\hat{\bf e} / (N-r({\bf X}))...
... & = & ( -13.6257 + (-.8391)(155.5175))/10, \\
& = & -14.4120.
\end{eqnarray*}


All of these values seem to be twice as large as the assumed values at the start. The REML process is iterative and the new estimates need to be used to re-form the MME, solve, and so forth.

These data were not very amenable to estimation of the maternal genetic or permanent environmental components because there were only two dams with progeny, and the estimation of the covariance between direct and maternal components would suffer because none of the animals appeared both as a dam and as an offspring. That is, the only link between direct and maternal effects is through relationships. To show that this is a problem, note that

\begin{displaymath}\hat{m}_{1} = -\frac{.4628}{3.24}(2.3295) = -.3328. \end{displaymath}

This relationship holds for nearly all animals in the analysis, except animals 3, 4, and 9. Thus, the estimate of the maternal genetic component for animal 1 is a function of the elements of ${\bf G}^{-1}$ and the solution for the direct genetic component. The correlation between $\hat{\bf a}$and $\hat{\bf m}$ was -.961, which is a lot stronger than the assumed prior correlation of -.196. There needs to be female animals with records that later have their own progeny in order to get sound REML estimates of the direct - maternal covariance. Thus, several generations of data are needed with good links between generations. To get a good estimate of the maternal ability of a bull's daughters, a bull has to have daughters which have had calves (i.e. grandprogeny). Often there are insufficient links between generations and this has given rise to many highly biased estimates of the direct - maternal genetic correlation.

Bayesian Estimation

The maternal effects model requires a slightly different approach because of the non-zero correlation between direct and maternal genetic effects. Let ${\bf g}' = ( {\bf a}' \ \ {\bf m}' )$, then

\begin{displaymath}{\bf g} \mid {\bf A},{\bf G} \sim N({\bf0}, ({\bf G}\otimes
{\bf A}) ), \end{displaymath}

and then

\begin{eqnarray*}p({\bf g} \mid {\bf A},{\bf G}) & = & \mid 2\pi ({\bf G}
\oti...
...ac{1}{2} {\bf g}'({\bf G}^{-1} \otimes {\bf A}^{-1})
{\bf g} ]
\end{eqnarray*}


Now look at

\begin{eqnarray*}{\bf g}'({\bf G}^{-1} \otimes {\bf A}^{-1})
{\bf g} & = & \lef...
...m} \end{array} \right), \\
& = & tr({\bf G}^{-1}{\bf S}_{g}),
\end{eqnarray*}


where

\begin{displaymath}{\bf S}_{g} = \left( \begin{array}{ll}
{\bf a}'{\bf A}^{-1}{...
...-1}{\bf a} & {\bf m}'{\bf A}^{-1}{\bf m}
\end{array} \right). \end{displaymath}

The prior distributions for the variance components are

\begin{eqnarray*}p(\sigma^{2}_{p} \mid v_{p},S^{2}_{p} ) & \propto &
(\sigma^{2...
... + 1)} \exp [ - \frac{1}{2}
tr({\bf G}^{-1}{\bf G}^{-1}_{o}) ],
\end{eqnarray*}


where vp, ve, $\nu$, S2p, S2e, and ${\bf G}_{o}$ are hyperparameters, p=2 is the dimension of ${\bf G}$. The prior distributions for $\sigma^{2}_{p}$ and $\sigma^{2}_{e}$are scaled inverted Chi-square distributions while that for ${\bf G}$is an inverse Wishart distribution, the multivariate equivalent of the inverse Chi-square distribution. Thus, we will need a good random Wishart matrix generator.

The joint posterior distribution is constructed and from this the fully conditional distributions for each component of the unknowns can be derived. The results are similar to the previous models. The vectors ${\bf a}$, ${\bf m}$, ${\bf p}$, and ${\bf b}$ in the MME are processed in the same manner as before with some small modifications. For the fixed effects,

\begin{eqnarray*}\hat{b}_{i} & = & ({\bf x}'_{i}{\bf x}_{i})^{-1}{\bf x}'_{i}
(...
...t{b}_{i} + RND*(\sigma^{2}_{e}/({\bf x}'_{i}
{\bf x}_{i}))^{.5}
\end{eqnarray*}


For the additive (direct) genetic effects,

\begin{eqnarray*}\hat{a}_{i} & = & ({\bf z}'_{1i}{\bf z}_{1i}+{\bf A}^{i,i}k_{11...
...2}_{e}/(
{\bf z}'_{1i}{\bf z}_{1i}+{\bf A}^{i,i}k_{11}))^{.5}.
\end{eqnarray*}


For the maternal genetic effects,

\begin{eqnarray*}\hat{m}_{i} & = & ({\bf z}'_{2i}{\bf z}_{2i}+{\bf A}^{i,i}k_{22...
...2}_{e}/(
{\bf z}'_{2i}{\bf z}_{2i}+{\bf A}^{i,i}k_{22}))^{.5}.
\end{eqnarray*}


For the maternal PE effects,

\begin{eqnarray*}\hat{p}_{i} & = & ({\bf z}'_{3i}{\bf z}_{3i}+k_{33})^{-1}
[ {\...
...RND*(\sigma^{2}_{e}/(
{\bf z}'_{3i}{\bf z}_{3i}+k_{33}))^{.5}.
\end{eqnarray*}


The computations are the same as before for $\sigma^{2}_{e}$ and $\sigma^{2}_{p}$, but are a little different for the genetic variances and covariance. For the residual component,

\begin{eqnarray*}{\bf e} & = & {\bf y}-{\bf Xb}-{\bf Z}_{1}{\bf a}-{\bf Z}_{2}{\...
...& = & ({\bf e}'{\bf e} + v_{e}S^{2}_{e}) /
CHI(\tilde{v}_{e}),
\end{eqnarray*}


for $\tilde{v}_{e} = (N+v_{e})$. For the maternal PE component,

\begin{displaymath}\sigma^{2}_{p} = ({\bf p}'{\bf p}+v_{p}S^{2}_{p}) /
CHI(\tilde{v}_{p}), \end{displaymath}

where $\tilde{v}_{p}=(n_{p}+v_{p})$, and np is the number of dams with offspring. For the genetic components, ${\bf G}$, the calculations are first,

\begin{displaymath}{\bf S}_{g} = \left( \begin{array}{ll}
{\bf a}'{\bf A}^{-1}...
...m}'{\bf A}^{-1}{\bf m}
\end{array} \right) + \nu {\bf G}_{o}, \end{displaymath}

then invert ${\bf S}_{g}$ and perform a Cholesky decomposition on the inverted matrix, i.e.,

\begin{displaymath}{\bf S}^{-1}_{g} = {\bf TT}', \end{displaymath}

then utilize a Wishart matrix generator such as

       CALL WISHRT(T,GI,ndf)
where ndf$=q+\nu$, then GI is the new ${\bf G}^{-1}$ to use in re-forming the MME for the next round of sampling. Thus,

\begin{displaymath}{\bf G} = ({\rm GI})^{-1}. \end{displaymath}

Then the next round of sampling begins. A critical point in this scheme is the Cholesky decomposition of ${\bf S}^{-1}_{g}$. A check should always be made that this calculation has been performed satisfactorily, i.e. that the routine did not encounter a 0 or negative number on the diagonals during the decomposition. However, the manner in which ${\bf S}_{g}$ is formed and presuming that ${\bf G}_{o}$ is positive definite, then this problem should not arise, but a check is still worthwhile.


next up previous

This LaTeX document is available as postscript or asAdobe PDF.

Larry Schaeffer
2000-03-17