This LaTeX document is available as postscript or asAdobe PDF.

Bayesian Estimation

These notes are based on Sorensen (1998) from a course on "Gibbs Sampling in Quantitative Genetics" given at AGBU, Armidale, NSW, Australia. The Gibbs sampling tool to derive estimates of variance components from the joint posterior distribution seems to be a very powerful tool to use with animal models. Computationally, one needs only a program that calculates solutions to Henderson's mixed model equations, very good random number generators, and computer time to process very large numbers of samples.

Begin with the simple single trait animal model. That is,

$${\bf y} = {\bf Xb}+{\bf Za} + {\bf e}.$$

The Bayesian approach is to derive the joint posterior distribution by application of Bayes theorem. If $\theta$ is a vector of random variables and ${\bf y}$ is the data vector, then

$$\begin{eqnarray*}p(\theta, {\bf y}) & = & p(\theta) \ p({\bf y} \mid \theta) \\ & = & p({\bf y}) \ p(\theta \mid {\bf y}) \end{eqnarray*}$$

Re-arranging gives

$$\begin{eqnarray*}p(\theta \mid {\bf y}) & = & \frac{p(\theta)p({\bf y} \mid \the... ... & = & {\rm posterior \ probability \ function \ of } \ \theta \end{eqnarray*}$$

In terms of the simple animal model, $\theta$ includes ${\bf b}$, ${\bf a}$, $\sigma^{2}_{a}$, and $\sigma^{2}_{e}$. The conditional distribution of ${\bf y}$ given $\theta$ is

$$\begin{eqnarray*}{\bf y} \mid {\bf b},{\bf a},\sigma^{2}_{a},\sigma^{2}_{e} & ... ...{\bf Za})' ({\bf y}-{\bf Xb}-{\bf Za})/2\sigma^{2}_{e} \right] \end{eqnarray*}$$

Prior distributions need to be assigned to the components in $\theta$, and these need to be multiplied together and times the conditional distribution of ${\bf y}$ given $\theta$. For the fixed effects vector, ${\bf b}$, there is little prior knowledge about the values that elements in that vector might have. This is represented by assuming

$$p({\bf b}) \propto {\rm constant}.$$

For ${\bf a}$, the vector of additive genetic values, quantitative genetics theory suggests that they follow a normal distribution, i.e.

$$\begin{eqnarray*}{\bf a} \mid {\bf A}, \sigma^{2}_{a} & \sim & N({\bf0}, {\bf ... ...left[ - {\bf a}'{\bf A}^{-1}{\bf a} / 2 \sigma^{2}_{a} \right], \end{eqnarray*}$$

where q is the length of ${\bf a}$. A natural estimator of $\sigma^{2}_{a}$ is ${\bf a}'{\bf A}^{-1}{\bf a} /q$, call it S²_a, which has

$$S^{2}_{a} \sim \chi^{2}_{q} \sigma^{2}_{a} /q.$$

Multiplying both sides by q and dividing by $\chi^{2}_{q}$ gives

$$\sigma^{2}_{a} \sim q S^{2}_{a} / \chi^{2}_{q}$$

which is a scaled, inverted Chi-square distribution, written as

$$p(\sigma^{2}_{a} \mid v_{a},S^{2}_{a}) \propto (\sigma^{2}_... ...t( -\frac{v_{a}}{2} \frac{S^{2}_{a}}{\sigma^{2}_{a}} \right),$$

where v_a and S²_a are hyperparameters with S²_a being a prior guess about the value of $\sigma^{2}_{a}$ and v_a being the degrees of belief in that prior value. Usually q is much larger than v_a and therefore, the data provide nearly all of the information about $\sigma^{2}_{a}$. Similarly, for the residual variance,

$$\begin{eqnarray*}p(\sigma^{2}_{e} \mid v_{e},S^{2}_{e}) & \propto & (\sigma^{2}... ... -\frac{v_{e}}{2} \frac{S^{2}_{e}}{\sigma^{2}_{e}} \right). \\ \end{eqnarray*}$$

Now form the joint posterior distribution as

$$p({\bf b},{\bf a},\sigma^{2}_{a},\sigma^{2}_{e} \mid {\bf y} ... ...({\bf y} \mid {\bf b},{\bf a},\sigma^{2}_{a}, \sigma^{2}_{e})$$

which can be written as

$$\propto (\sigma^{2}_{e})^{-( \frac{N+v_{e}}{2}+1)} \exp \l... ...bf Za})'({\bf y}-{\bf Xb}-{\bf Za}) + v_{e}S^{2}_{e}) \right]$$

$$(\sigma^{2}_{a})^{-(\frac{q+v_{a}}{2}+1)} \exp \left[ - \fra... ...{a}} ({\bf a}'{\bf A}^{-1}{\bf a} + v_{a}S^{2}_{a}) \right].$$

In order to implement Gibbs sampling, all of the fully conditional posterior distributions (one for each component of $\theta$ ) need to be derived from the above joint posterior distribution. Let

$$\begin{eqnarray*}{\bf W} & = & ( {\bf X} \ \ {\bf Z} ), \\ \beta' & = & ( {\bf ... ...W}'{\bf W} + \Sigma \\ {\bf C}\hat{\beta} & = & {\bf W}'{\bf y} \end{eqnarray*}$$

Now we adopt a new notation, let

$$\beta' = ( \beta_{i} \ \ \beta'_{-i} ),$$

where $\beta_{i}$ is a scalar representing just one element of the vector $\beta$, and $\beta_{-i}$ is a vector representing all of the other elements except $\beta_{i}$. Similarly, ${\bf C}$ and ${\bf W}$ can be partitioned in the same manner as

$$\begin{eqnarray*}{\bf W}' & = & ( {\bf W}_{i} \ \ {\bf W}_{-i} )' \\ {\bf C} & ... ...{i,-i} \\ {\bf C}_{-i,i} & {\bf C}_{-i,-i} \end{array} \right). \end{eqnarray*}$$

In general terms, the conditional posterior distribution of $\beta$ is

$$\begin{eqnarray*}\beta_{i} \mid \beta_{-i}, \sigma^{2}_{a}, \sigma^{2}_{e},{\bf y} & \sim & N(\hat{\beta}_{i}, C^{-1}_{i,i} \sigma^{2}_{e}) \end{eqnarray*}$$

where

$$C_{i,i}\hat{\beta}_{i} = ({\bf W}'_{i}{\bf y} - {\bf C}_{i,-i} \beta_{-i}).$$

Then

$$b_{i} \mid {\bf b}_{-i}, {\bf a}, \sigma^{2}_{a}, \sigma^{2}_... ...{\bf y} \ \sim \ N(\hat{b}_{i}, C^{-1}_{i,i} \sigma^{2}_{e} ),$$

for

$$C_{i,i} = {\bf x}'_{i}{\bf x}_{i}.$$

Also,

$$a_{i} \mid {\bf b},{\bf a}_{-i}, \sigma^{2}_{a}, \sigma^{2}_{... ... {\bf y} \ \sim \ N(\hat{a}_{i}, C^{-1}_{i,i} \sigma^{2}_{e}),$$

where

$$C_{i,i} = ( {\bf z}'_{i}{\bf z}_{i} + A^{i,i}k),$$

for $k=\sigma^{2}_{e}/ \sigma^{2}_{a}$. The conditional posterior distributions for the variances are

$$\sigma^{2}_{a} \mid {\bf b},{\bf a}, \sigma^{2}_{e},{\bf y} \sim \tilde{v}_{a} \tilde{S}^{2}_{a} \chi^{-2}_{\tilde{v}_{a}}$$

for $\tilde{v}_{a} = q + v_{a}$, and $\tilde{S}^{2}_{a} = ({\bf a}'{\bf A}^{-1}{\bf a} + v_{a}S^{2}_{a} ) / \tilde{v}_{a},$and

$$\sigma^{2}_{e} \mid {\bf b},{\bf a}, \sigma^{2}_{a},{\bf y} \sim \tilde{v}_{e} \tilde{S}^{2}_{e} \chi^{-2}_{\tilde{v}_{e}}$$

for $\tilde{v}_{e} = N + v_{e}$, and $\tilde{S}^{2}_{e} = ({\bf e}'{\bf e} + v_{e}S^{2}_{e} ) / \tilde{v}_{e}$, and ${\bf e} = {\bf y}-{\bf Xb}-{\bf Za}$.

Computational Scheme

Gibbs sampling is much like Gauss-Seidel iteration, except that

1.

when a new solution is calculated a random amount is added to it based upon its conditional posterior distribution variance before proceeding to the next equation, and

2.

after all equations have been processed, new values of the variances are calculated and a new variance ratio is determined prior to beginning the next round.

Suppose we have the following MME for five animals:

$$\left( \begin{array}{rrrrrr} 5 & 1 & 1 & 1 & 1 & 1 \\ 1 & ... ...2 \\ 38.5 \\ 48.9 \\ 64.3 \\ 50.5 \\ 36.0 \end{array} \right),$$

where $k=\sigma^{2}_{e}/\sigma^{2}_{a} = 14$, and

$${\bf A}^{-1} = \frac{1}{14} \left( \begin{array}{rrrrr} 28 &... ...14 & 36 & -16 \\ 0 & -16 & 0 & -16 & 32 \end{array} \right).$$

Let the starting values for $\beta = \left( \begin{array}{rrrrrr} 0 & 0 & 0 & 0 & 0 & 0 \end{array} \right)$, and let v_a=v_e=10, and $S^{2}_{e}=93\frac{1}{3}$ and $S^{2}_{a}=6\frac{2}{3}$, so that k=14. Let RND represent a random normal deviate from a random normal deviate generator, and let CHI(idf) represent a random Chi-square variate from a random Chi-Square variate generator. To begin, let $\sigma^{2}_{e}=S^{2}_{e}$and $\sigma^{2}_{a}=S^{2}_{a}$. Below are descriptions of calculations in the first two rounds.

Round 1

$$\begin{eqnarray*}\hat{\mu} & = & (238.2 - a_{1}-a_{2}-a_{3}-a_{4}-a_{5})/5 \\ ... ...{.5} \\ & = & -.9316 + (-.6472)(1.6817) \\ & = & -2.0200 \\ \end{eqnarray*}$$

Now calculate the residuals and their sum of squares in order to obtain a new residual variance.

$$\begin{eqnarray*}e_{1} & = & 38.5 - 42.41 -1.9067 = -5.8167 \\ e_{2} & = & 48.9... ... 42.41 + 2.0200 = -4.3900 \\ {\bf e}'{\bf e} & = & 705.1503 \\ \end{eqnarray*}$$

The new residual variance is

$$\begin{eqnarray*}\sigma^{2}_{e} & = & ({\bf e}'{\bf e} + v_{e}S^{2}_{e}) / CHI(1... ...\ & = & (705.1503 + (10)(93.3333))/17.1321 \\ & = & 95.6382. \end{eqnarray*}$$

The additive genetic variance requires calculation of ${\bf a}'{\bf A}^{-1}{\bf a}$ using the a-values obtained above, which gives

$${\bf a}'{\bf A}^{-1}{\bf a} = 19.85586.$$

Then

$$\begin{eqnarray*}\sigma^{2}_{a} & = & ( {\bf a}'{\bf A}^{-1}{\bf a} + v_{a} S^{... ...\ & = & ( 19.85586 + (10)(6.66667))/10.7341 \\ & = & 8.0605. \end{eqnarray*}$$

The new variance ratio becomes

k = 95.6382 / 8.0605 = 11.8650.

Round 2

Round 2 begins by re-forming the MME using the new variance ratio. The equations have changed to

$$\left( \begin{array}{rrrrrr} 5 & 1 & 1 & 1 & 1 & 1 \\ 1 & ... ...2 \\ 38.5 \\ 48.9 \\ 64.3 \\ 50.5 \\ 36.0 \end{array} \right).$$

Now the process is repeated using the last values of $\mu$ and ${\bf a}$ and $\sigma^{2}_{e}$.

$$\begin{eqnarray*}\hat{\mu} & = & (238.2 - a_{1}-a_{2}-a_{3}-a_{4}-a_{5})/5 \\ ... ...\\ a_{5} & = & -2.6253 + (.8184)(1.8442) \\ & = & -1.1160 \\ \end{eqnarray*}$$

The residuals and their sum of squares are

$$\begin{eqnarray*}e_{1} & = & 38.5 - 51.41 +1.3999 = -11.5101 \\ e_{2} & = & 48.... ...51.41 + 1.1160 = -14.2940 \\ {\bf e}'{\bf e} & = & 627.1630 \\ \end{eqnarray*}$$

The new residual variance is

$$\begin{eqnarray*}\sigma^{2}_{e} & = & ({\bf e}'{\bf e} + v_{e}S^{2}_{e}) / CHI(1... ...\ & = & (627.1630 + (10)(93.3333))/20.4957 \\ & = & 76.1377. \end{eqnarray*}$$

The additive genetic variance is

$$\begin{eqnarray*}\sigma^{2}_{a} & = & ( {\bf a}'{\bf A}^{-1}{\bf a} + v_{a} S^{... ...\\ & = & ( 36.8306 + (10)(6.66667))/16.6012 \\ & = & 6.2343. \end{eqnarray*}$$

The new variance ratio becomes

k = 76.1377 / 6.2343 = 12.2127.

Continue taking samples for thousands of rounds.

Burn In Period

Generally, the Gibbs chain of samples does not immediately converge to give samples from the joint posterior distribution. A burn in period is usually allowed during which the sampling process moves from the initial values of the parameters to those from the joint posterior distribution. The length of the burn-in period is usually judged by visually inspecting a plot of sample values across rounds. In the example figure below, the burn-in period should likely be 2000 rounds. This figure represents an animal model analysis of over 323,000 records on slightly over 1 million animals. The trait was one of the type traits in the Holstein breed. The samples generated during the burn-in period are typically discarded from further use.

Estimates

Each round of Gibbs sampling is dependent on the results of the previous round. Depending on the total number of observations and parameters, one round may be positively correlated with the next twenty to three hundred rounds. Even though samples are correlated, the user can determine the effective number of samples in the total by calculating lag correlations. Suppose a total of 12,000 samples (after removing the burn-in rounds) gave an effective number of samples equal to 500. This implies that every 24^th sample should be uncorrelated.

An overall estimate of a parameter can be obtained by averaging all of the 12,000 samples (after the burn-in). However, to derive a confidence interval or to plot the distribution of the samples or to calculate the standard deviation of the estimate, the variance of the 500 independent samples should be used.

Influence of the Priors

In the small example, v_a=v_e=10 whereas N was only 5. Thus, the prior values of the variances received more weight than information coming from the data. This is probably appropriate for this small example, but if N were 5,000,000, then the influence of the priors would be next to none. The amount of influence of the priors is not directly determined by the ratio of v_i to N. In the small example, even though $v_{a}/(N+v_{a})=\frac{2}{3}$, the influence of S²_a could be greater than $\frac{2}{3}$. Schenkel (1998) applied some methods recommended to him by Gianola in his thesis work.

Long Chain or Many Chains?

Early papers on MCMC (Monte Carlo Markov Chain) methods recommended running many chains of samples and then averaging the final values from each chain. This was to insure independence of the samples. That philosophy seems to have changed to a recommendation of one single long chain. For animal breeding applications this could mean 100,000 samples or more. If a month is needed to run 50,000 samples, then maybe three chains of 50,000 would be preferable. If only an hour is needed for 50,000 samples, then 1,000,000 samples would not be difficult to run. The main conclusion is that it is very difficult to know if or when a Gibbs chain has converged to the joint posterior distribution. This is currently an active area of statistical research.

This LaTeX document is available as postscript or asAdobe PDF.