1. No category

Technical note: How to use Winbugs to draw inferences in animal models
L. H. Damgaard
J ANIM SCI 2007, 85:1363-1368.
doi: 10.2527/jas.2006-543 originally published online January 3, 2007
The online version of this article, along with updated information and services, is located on
the World Wide Web at:
http://www.journalofanimalscience.org/content/85/6/1363
www.asas.org
Downloaded from www.journalofanimalscience.org by guest on September 19, 2014
Technical note: How to use Winbugs to draw inferences in animal models1
L. H. Damgaard*†2
*Department of Genetics and Biotechnology, Danish Institute of Agricultural Sciences, PO Box 50, 8830 Tjele,
Denmark; and †Department of Natural Sciences, Royal Veterinary and Agricultural University,
Thorvaldsensvej 40, 1870 Frederiksberg C, Denmark
ABSTRACT: This paper deals with Bayesian inferences of animal models using Gibbs sampling. First,
we suggest a general and efficient method for updating
additive genetic effects, in which the computational cost
is independent of the pedigree depth and increases linearly only with the size of the pedigree. Second, we
show how this approach can be used to draw inferences
from a wide range of animal models using the computer
package Winbugs. Finally, we illustrate the approach
in a simulation study, in which the data are generated
and analyzed using Winbugs according to a linear
model with i.i.d errors having Student’s t distributions.
In conclusion, Winbugs can be used to make inferences
in small-sized, quantitative, genetic data sets applying
a wide range of animal models that are not yet standard
in the animal breeding literature.
Key words: Bayesian inference, Winbugs, genetics
©2007 American Society of Animal Science. All rights reserved.
INTRODUCTION
The beginning point for the development of new models in quantitative genetics often originates with a fixed
model that assumes individual animals to be independent. The idea is to include the effect of additive gene
action on a linear additive scale, together with other
explanatory variables, and then to assume a multivariate normal distribution for the unobserved additive genetic effects. Examples include the mixed model for
ordered categorical data (Sorensen et al., 1995), which
is based on the threshold liability concept first proposed
by Wright (1934), and the mixed model for survival
data (Ducrocq and Casella, 1996), which is based on
the Cox proportional hazards model (Cox, 1972).
Although well-documented software often exists for
the fixed analog of the additive genetic model, and in
some cases also a simple mixed model, the applicability
of such software in animal breeding is often limited.
The reason is that it is impossible to set up the additive
genetic relationship matrix for models with a pedigree
structure more complex than that of a sire model.
Hence, the limiting factors for applying new models
suggested in the statistical literature to analyze quanti-
J. Anim. Sci. 2007. 85:1363–1368
doi:10.2527/jas.2006-543
tative genetic data are the time and resources that necessarily must be allocated to the development of new
software.
In this paper the goal was first to develop a computationally efficient method for updating additive genetic
effects, and second to show how this method could be
used to analyze genetic data using Bayesian probability
models with the flexible computer package Winbugs
(Spiegelhalter et al., 2006). The idea was first proposed
by Mrode and Thompson (1989), in which they parameterized the mixed linear model in terms of the transformed additive genetic effects based on the usual partition of the additive genetic relationship matrix (Henderson, 1976). Finally, the approach was illustrated in
a simulation study, in which the data were generated
according to a linear model, with i.i.d errors having
Student’s t distributions.
MATERIALS AND METHODS
Animal Care and Use Committee approval was not
obtained for this study because only no animals were
used.
The Bayesian Additive Genetic Model
1
I thank Ignacy Misztal for useful discussions on the partition of
the additive genetic relationship matrix and associated calculations.
2
Corresponding author: [email protected]
Received August 8, 2006.
Accepted December 21, 2006.
Bayesian probability models are constructed by specifying the joint distribution, say p(y,θ), of the observable
y, which is the data, and the unobservable θ, which
comprises the model parameters, missing data, and so
on. As soon as the data are observed, model inferences
1363
Downloaded from www.journalofanimalscience.org by guest on September 19, 2014
1364
Damgaard
are based on the posterior distribution p(y | θ). Using
Bayes’s theorem, the posterior distribution is, up to
proportionality, given by the product of the sampling
distribution (likelihood) p(y | θ) and the prior p(θ):
p(θ | y) ∝ p(y | θ)p(θ).
For the underlying genetic model, we assume an additive genetic infinitesimal model (Bulmer, 1980). We will
classify the parameter vector into a vector of additive
genetic effects, a, and associated covariance matrix, G,
and the remaining parameters, η. Moreover, we adapt
the often-used assumption of conditional independence,
which implies that the data are independent given the
model parameters. A priori, we assume that p(θ) =
p(η)p(a | G)p(G). According to quantitative genetic theory, a | G is assumed to be multivariate normally distributed, with a mean vector 0 and a covariance matrix
G ⊗ A, where A is the additive genetic relationship
matrix. The prior specification for η and G will depend
on the problem at hand and will be left unspecified until,
in the last part of this paper, an example is considered.
Under the above assumption, a Bayesian additive genetic model may be represented by
p(θ | y)
∝
[1]
Π p(yi | η, a)⎤⎥⎦p(η)Na|G (0, G ⊗ A)p(G).
⎡
⎢
⎣ i
Reparameterized Model
In the following model, the additive genetic effects
are reparameterized according to a = (L ⊗ TD1/2)γ,
where L is the lower triangular Cholesky decomposition
of G (G = LL′), and T and D correspond to the factorization A = TDT′ (Henderson, 1976). This implies that the
vector γ is standard, multivariate normally distributed.
By the rules for the multidimensional change of random
variables, the posterior distribution for the reparameterized model is, up to proportionality, given by
p(θ˜ | y)
∝
[2]
Π p(yi | η, a) | a=(L⊗TD1/2)γ⎤⎥⎦p(η)Nγ (0 I)p(L),
⎡
⎢
⎣ i
where a has been exchanged by (L ⊗ TD1/2)γ, and the
prior distribution for L is given by the prior distribution
for G and the transformation G = LL′. For G to be
positive definite, the diagonal elements of L must be
greater than zero.
The parameter vector θ˜ is given by (η,γ,L). Therefore,
if Markov chain Monte Carlo methods are used for the
inferences, then posterior samples are obtained of γ
and L (Metropolis et al., 1953; Hastings, 1970). These
are, however, easily transformed to posterior samples
of a and G, and Bayesian inferences including posterior
point estimates and credibility intervals are done as
usual.
It is important to note that the prior distribution is
not generally invariant to the reparameterization. As
an example, consider the situation in which we have
only a direct additive genetic effect, set G = σa2, from
which it follows that L = √σa2. If here we assume a priori
a bounded uniform prior for L on the interval (a, b)
with 0 < a < b, then the prior distribution for σa2 is given
1
by (σa2)−1/2/(b2 − a2) on the interval (a2,b2), which is
2
bounded but not uniform.
Gibbs Updating
In a Gibbs sampler, posterior samples are obtained
by repeatedly updating from the full conditionals
(Sorensen and Gianola, 2002). These are easily derived,
at least up to proportionality, by retaining the parameters of interest from the posterior distribution.
From the reparameterized posterior distribution
(equation [2]), it follows that updating from the full
conditional for γ requires the computation of (L ⊗
TD1/2)γ. Because the element Tij of T is the fraction of
genes of animal i inherited by animal j, the average
number of nonzero elements in one row of T is the
average number of ancestors per animal. The T matrix
is therefore relatively sparse, and it would be an inefficient approach to form T and evaluate (L ⊗ TD1/2)γ
explicitly. Fortunately, Quaas (1989) presented a recursive algorithm for calculating v = Ts (where s is an
arbitrary vector), working directly from a list of sires
and dams and therefore with cost independent of the
pedigree depth. Thus, the matrix T is never formed.
Assuming that animals have been sorted from the oldest to the youngest and renumbered from 1 to q (parents
precede offspring), the recursive formula is
vi = si + (vsi + vdi)/2
for i = 1, ..., q,
[3]
where the subscripts si and di are the sire and dam
index of animal i. If a sire (dam) is unknown, then vsi =
0 (vdi = 0). This result easily extends to the calculation
of a = (L ⊗ TD1/2)γ.
Clearly, this result is useful for updating the elements of γ = (γ1, ..., γq), irrespective of whether it is
done jointly or univariately. We will focus on the latter,
because this allows us to use the standard computer
package Winbugs. In Table 1, the general idea of how
to update the elements γ1, ..., γq for the univariate case
is sketched. A specific example of how to fit an animal
model in Winbugs is given in the next section. For ease
of representation, one trait with a direct additive genetic effect is considered, in which case L is a scalar
and equal to √σa2. Extension to several additive genetic
effects and multivariate models is straightforward but
notationally cumbersome and will not be considered
here.
Downloaded from www.journalofanimalscience.org by guest on September 19, 2014
1365
Gibbs updating of additive genetic effects
Table 1. General subroutine for one round of univariate
updating of transformed additive genetic effects1
INITIALIZE working vector v
DO i = 1, q
IF (animal i has no record)
UPDATE γi FROM
portant because computational advances will enable
the fitting of increasingly complex models.
EXAMPLE
Student’s t Animal Model
⎛ 1 ⎞
γi | rest ∝ exp⎜− γi2⎟
⎝ 2 ⎠
1
CALCULATE νi = √Dii × γi (νsi + νdi)
2
ELSE IF (animal i has record)
UPDATE γi FROM
γi | rest ∝ p(yg(i) | η, ai)|a
i
=
⎡
√ ⎢⎣√
2
a
σ
⎤
1
D × γ + (ν + ν ) ⎥
a
i 2 si
di ⎦
⎛ 1 ⎞
exp⎜− γi2 ⎟
⎝ 2 ⎠
1
CALCULATE νi = √Dii × γi (νsi + νdi)
2
END IF
END DO
1
γi = the transformed genetic effect; y = the observation; D = the
diagonal element of the D matrix; si = the sire id; di = the dam id;
v = the working vector; g(i) = the index function that relates the
genetic effect to the observation; q = the number of animals in the
pedigree.
Different strategies for univariate updating of parameters can be applied if the normalizing constant
cannot be found. These include adaptive rejection sampling (Gilks and Wild, 1992; Gilks and Best, 1995) and a
Metropolis-Hastings algorithm (Metropolis et al., 1953;
Hastings, 1970). In situations in which it is possible to
find the normalizing constant, direct updating using
standard methods is an obvious choice. Winbugs frees
the user from these computational details. This advantage makes Bayesian inference applicable for a wider
audience. However, it is still the analyst’s responsibility
to ensure that the model assumptions are satisfied;
otherwise, misleading inferences may be obtained. In
the future, model validation will be even more im-
To illustrate the method, we generate data according
to an animal model with errors independently and identically distributed as t(0, σ2, ν). Here, σ2 is a positive
scale parameter and ν is a degrees of freedom parameter. If ν is ≤2, then the variance of the Student’s t pro⎛ ν
⎞
σ2⎟ is infinite. When ν goes toward infinity,
cess ⎜
⎝ν − 2
⎠
the Student’s t process tends toward a normal process
(e.g., t(0, σ2, ν) → N(0, σ2) for ν → ∞). Compared with
the normal distribution, the Student’s t distribution has
thicker tails and is therefore less sensitive to outliers,
which may have a marked impact on inferences in the
normal model (Zellner, 1976). The use of t distributed
errors in models for analyzing quantitative genetic data
was addressed only recently (Strande´n and Gianola
1999). I am not aware of any easy-to-use software for
inferring this model, and the possibility of using Winbugs is therefore relevant.
Let Y = (Y1,...,Yn) denote the random vector of observations. The sampling distribution is given by
n
p(y | θ) =
Π t(xi′ β + zi′ a, σ2, ν),
where xi′ is the ith row of the n × p design matrix X, β
is the p × 1 vector of fixed effects, zi′ is the ith row of
the n × q design matrix Z, and a is the q × 1 vector of
additive genetic effects that, according to quantitative
genetic theory, is assumed to be multivariate normally
distributed with mean vector 0 and covariance matrix
Table 2. Data specification in Winbugs1
Data file 1
X[]
0
0
—
1
2
—
1
END
Data file 2
nd[]
1
2
—
200
Y[]
0
0
—
100.8
91.4
—
106.6
Sq_D[]
1
1
—
0.75
0.50
—
0.45
[4]
i=1
ID[]
1
2
—
201
202
—
2,284
SID[]
0
0
—
1
1
—
1,856
DID[]
0
0
—
0
3
—
1,911
ad[]
201
202
—
—
2,284
END
1
Data file 1: X = the covariate; Y = the observation; sq_D = the square root of diagonal element of the D
matrix; SID = the sire id; DID = the dam id. Data file 2: nd = the pointer to observation number for an
animal without a record; ad = the pointer to observation number for an animal with a record.
Downloaded from www.journalofanimalscience.org by guest on September 19, 2014
1366
Damgaard
Table 3. Model specification in Winbugs for a Student’s t animal model1
# Specification of the reparameterized sampling distribution
1. model {
2. for( i in 1 : 2084 ) {
3.
Y[ ad[i] ] ∼ dt( mu[ ad[i] ], tau, ndf )
4.
mu[ ad[i] ] ← beta[ X[ ad[i] ] ] + v[ ad[i] ] * gam.tau
5.
v[ ad[i] ] ← gam[ ID[ ad[i] ] ] * sq_D[ ad[i] ] + ( v[ SID[ ad[i] ] ] + v[ DID[ ad[i] ] ] )/2.0
6. }
7. for ( i in 1 : 200 ){
8.
v[ nd[i] ] ← gam[ ID[ nd[i] ] ] * sq_D[ nd[i] ] + ( v[ SID[ nd[i] ] ] + v[ DID[ nd[i] ] ] )/2.0
9. }
# Working elements (zero) used for animals without one or both parents
10.
v[2285] ← 0.0
11.
v[2286] ← 0.0
# Specification of prior distributions
12.
for (i in 1 : 2284 ){
13.
gam[i] ∼ dnorm( 0, 1.0 )
14.
}
15.
ndf∼dunif( 2, 100 )
16.
beta[1] ∼ dnorm( 0.0, 0.001 )
17.
beta[2] ∼ dnorm( 0.0, 0.001 )
18.
tau ∼ dgamma( 0.001, 0.001 )
19.
gam.tau ∼ dgamma( 0.001, 0.001 )
# Specification of functions of model parameters of inferential interest
20.
var_a ← gam.tau * gam.tau
21.
var_e ← 1 / ( tau )
22. }
1
Y = the observation; X = the covariate; sq_D = the square root of the diagonal element of the D matrix;
SID = the sire id; DID = the dam id; v = the working vector; nd = the pointer to observation number for an
animal without a record; ad = the pointer to observation number for an animal with a record; gam = the
transformed additive genetic effect; ndf = the degrees of freedom parameter; beta = the regression parameter;
tau = the reciprocal scale parameter; gam.tau = the square root of genetic variance; var_a = the genetic
variance; var_e = the scale parameter. Winbugs notation: dt = the t-density; dnorm = the normal density;
dunif = the uniform density; and dgamma = the gamma density.
Aσa2. The matrix A is the additive genetic relationship
matrix and σa2 is the additive genetic variance in the
base population. θ is the complete parameter vector,
and a priori we assume that β1, ..., βp, (a, σa2), σ2, and
ν are independent.
Data Simulation
The simulation study was designed to mimic a pedigree structure from a potential study in pig breeding.
More specifically, the pedigree was generated by selecting a multiplier herd from Danbred’s database, which
is the major pig breeding organization in Denmark, and
tracing the sows with first farrowing in the period from
January 1, 2002, to December 31, 2002, back 3 generations. Altogether, the pedigree included 2,284 animals
and the records were simulated for the 2,084 oldest
animals according to model [4]. In addition to σ2 = 100
and ν = 10, the model included a fixed effect with 2
levels equally assigned to the animals (β1 = 100, β2 =
120), and an animal effect with additive genetic variance (σa2 = 30).
Data Analysis
The data were analyzed with the same model as used
to generate the data. The assumptions were vague, nor-
mally distributed priors for fixed effects [N(0,1000)],
gamma distributed priors for the additive genetic variance
component
and
the
scale
parameter
[gamma(0.001,0.001)], and a bounded uniform prior for
the degrees of freedom parameter [uniform(2,100)].
By studying the instructive examples given in the
manual for Winbugs (Spiegelhalter et al., 2006), it follows that data and model specification can be done in
various ways. Data for the present analysis are given
in 2 files (Table 2). The first data file included the actual
data for animals without and with records, in which
a zero indicated missing data or an unknown parent.
Because the “if-then” statement does not exist in Winbugs and animals without and with records are handled
differently, the second data file included 2 vectors that
pointed to the position of the animals without and with
records in the first data file.
Except for the line numbers, Table 3 is a copy of the
model specification in Winbugs that was used to analyze
the simulated data. The model was fitted by running 3
independent chains for a total of 10,000 iterations, in
which the model parameters were saved every 10th
iteration. After the initial values for levels of fixed effects, genetic variance, scale parameter, and degree of
freedom parameter were assigned, initial values for additive genetic effects were generated by Winbugs. Different beginning values were used in each chain (see
Downloaded from www.journalofanimalscience.org by guest on September 19, 2014
1367
Gibbs updating of additive genetic effects
Table 4. Marginal posterior summary statistics1
Variable
Chain
Begin. val.
True
Mean
SD
95% HPD
β1
β1
β1
C1
C2
C3
90
120
80
100
—
—
99.8
99.8
19.8
0.60
0.60
0.61
98.1–101.8
98.6–101.1
98.7–101.1
β2
β2
β2
C1
C2
C3
90
80
90
120
—
—
119.9
119.9
119.9
0.60
0.60
0.60
118.0–120.9
118.7–121.1
118.7–121.1
σa2
C1
5
30
31.2
6.1
15.7–42.4
σa2
C2
50
—
31.2
6.0
20.1–43.9
σa2
C3
20
—
31.5
6.2
20.0–44.1
σ2
σ2
σ2
C1
C2
C3
200
100
50
100
—
—
101.9
102.7
100.9
7.5
7.7
7.7
73.0–115.4
88.4–117.4
85.7–116.4
ν
ν
ν
C1
C2
C3
5
80
30
10
—
—
13.0
13.4
12.6
4.3
5.3
4.7
5.4–24.7
6.8–27.4
6.9–26.9
1
β = the regression parameter; σa2 = the additive genetic variance component; σ2 = the scale parameter;
ν = the degrees of freedom parameter; C = the Markov chain Monte Carlo chain, Begin. val. = the beginning
values; 95% HPD = the 95% highest posterior density region.
Table 4). Although convergence for additive genetic effects was not checked, it was concluded, based on plots
of the Gibbs chains for the remaining parameters, that
an acceptable degree of convergence had been obtained.
This was further supported by the fact that posterior
summary statistics were very similar for the 3 chains
(Table 4). The scale parameter and the degree of freedom parameter were the most poorly behaved parameters in terms of mixing, with an autocorrelation at lag
20 (saved iterations) equal to 0.18 and 0.44, respectively, for all 3 chains. Note that those parameters also
are the parameters for which the posterior summary
statistics varied the most, and one may consider extending the Gibbs chain to decrease the Monte Carlo
error. For posterior summarization, the first 100 saved
iterations were discarded and results for the 3 chains
are given in Table 4.
Model Extensions
A clear advantage with the model specification in
Winbugs is that very few modifications are necessary
if one is interested in fitting alternative models. If, for
example, we want to run an animal model with normal
errors, then we only need to change “dt(.)” in line 3 to
“dnorm(.)” in line 3. If we want to fit additional fixed
or random effects, then these should be included on line
4, and an obvious place to include their prior specification is after line 17. Finally, if we want to estimate
heritability, then the formula for this variable can be
defined after line 22. Note that the second argument
in the “dnorm” statement, which is the code for a normal
distribution, is not the variance, but is one over the
variance.
In the post-Gibbs analysis, one may take advantage
of the tools available in Winbugs. Alternatively, one
can use the facility for outputting Gibbs samples in a
format that is compatible with the post-Gibbs analysis
computer package, CODA (Best et al., 1997).
Results
The posterior mean values and 95% highest posterior
density regions given in Table 4 show that the parameters can be correctly inferred using the proposed approach in Winbugs. In all cases, the 95% highest posterior density regions assign relatively high probability
mass to values of the parameters in the neighborhood
of the true values.
DISCUSSION
A general and computationally efficient method for
updating additive genetic effects in Bayesian inferences
of animal models using Gibbs sampling was presented.
It also was shown how the method can be used to draw
inferences from animal models using the flexible computer package Winbugs (Spiegelhalter et al., 2006).
The idea sketched here for updating additive genetic
effects is computationally efficient and relevant not only
for fitting models in Winbugs. The reason is that the
cost associated with updating additive genetic effects
is independent of the pedigree depth and increases linearly only with the size of the pedigree. The idea of
parameterizing an animal model in terms of transformed additive genetic effects is not new (Mrode and
Thompson, 1989). In a Bayesian implementation of a
genetically structured residual variance model, Sorensen and Waagepetersen (2003) used Langevin-Hastings
(Besag, 1994) to jointly update transformed additive
genetic effects. According to these authors, it is easier
to obtain a well-mixing Markov chain Monte Carlo algo-
Downloaded from www.journalofanimalscience.org by guest on September 19, 2014
1368
Damgaard
rithm for the posterior distribution of the transformed
variables than the original genetic effects. We have had
a similar experience in a Bayesian implementation of
the semiparametric Cox model with time varying genetic effects (Damgaard, 2006).
The flexibility of Winbugs is best illustrated by looking at the long list of examples provided with the documentation. Among many others, these examples include models for continuous data, counting data, binary
and categorical data, and survival data. Following the
example in this paper, it should be straightforward to
extend these models to animal models and use Winbugs
to make additive genetic inferences. In addition to the
possibility of fitting a wide range of standard Bayesian
models, it is also possible for the user to specify his or
her own distribution at each step of Bayesian model
building.
As soon as data have been set up for Winbugs, it is
very easy to fit alternative models. This feature enables
the analyst to consider a large group of competing models and use Bayesian model selection criteria to select
the best one for drawing inferences (e.g., Sorensen and
Gianola, 2002). This is clearly an advantage in that it
will decrease the risk of using an incorrect or inadequate model to interpret data. The reason is that a
larger group of models will tend to be considered, compared with a situation in which resources must be allocated to implement new software to fit alternative
models.
Unfortunately, the use of Winbugs in quantitative
genetics is limited to relatively small data sets, beginning from a few thousand records. The analysis performed in this paper took about 1 h on a Pentium computer (1.6 GHz, 512 MB RAM). Hence, it is unlikely
that we will see Winbugs used for genetic evaluation. Its
applicability is limited to inference in small populations
and experiments designed with the number of records,
for example, limited by a high recording cost, such as
in many behavioral studies. In addition, Winbugs may
in turn be useful for educational purposes.
LITERATURE CITED
Besag, J. 1994. Contribution to the discussion paper by Grenander
and Miller. J. R. Statist. Soc. B. 56:591–592.
Best, N. G., M. K. Cowles, and S. K. Vines. 1997. CODA: Convergence
diagnosis and output analysis software for Gibbs sampling output. Available: http://www.mrc-bsu.cam.ac.uk/bugs/classic/
coda04/readme.shtml. Accessed Aug. 5, 2006.
Bulmer, M. G. 1980. The Mathematical Theory of Quantitative Genetics. Clarendon Press, Oxford, UK.
Cox, D. R. 1972. Regression models and life tables. J. R. Statist. Soc.
B 34:187–220.
Damgaard, L. H. 2006. Joint quantitative genetic analysis of survival,
linear Gaussian and ordered categorical traits. Proc. 8th World
Congr. Genet. Appl. Livest. Prod., Belo Horizonte, Brazil. CDROM communication No. 2606.
Ducrocq, V., and G. Casella. 1996. A Bayesian analysis of mixed
survival models. Genet. Sel. Evol. 28:505–529.
Gilks, W., and R. N. G. Best. 1995. Adaptive rejection Metropolis
sampling within Gibbs sampling. Appl. Stat. 44:455–472.
Gilks, W., and R. P. Wild. 1992. Adaptive rejective sampling for Gibbs
sampling. Appl. Stat. 41:337–348.
Hastings, W. K. 1970. Monte Carlo sampling methods using Markov
chains and their application. Biometrika 57:97–109.
Henderson, C. R. 1976. A simple method for the inverse of a numerator
relationship matrix used in prediction of breeding values. Biometrics 32:69–83.
Metropolis, N. A., W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and
E. Teller. 1953. Equations of state calculations by fast computing
machines. J. Chem. Phys. 21:1087–1092.
Mrode, R., and R. Thompson. 1989. An alternative algorithm for
incorporating the relationships between animals in estimating
variance components. J. Anim. Breed. Genet. 106:89–95.
Quaas, R. L. 1989. Transformed mixed model equations: A recursive
algorithm to eliminate A−1. J. Dairy Sci. 72:1937–1941.
Sorensen, D., S. Andersen, D. Gianola, and I. R. Korsgaard. 1995.
Bayesian inference in threshold models using Gibbs sampling.
Genet. Sel. Evol. 27:229–249.
Sorensen, D., and D. Gianola. 2002. Likelihood, Bayesian, and MCMC
Methods in Quantitative Genetics. Springer, New York, NY.
Sorensen, D., and R. Waagepetersen. 2003. Normal linear models
with genetically structured residual variance heterogeneity: A
case study. Genet. Res. Camb. 82:207–222.
Spiegelhalter, D., A. Thomas, N. Best, and D. Lunn. 2006. Winbugs
1.4.1. Bayesian Inference Using Gibbs Sampling. Available:
http://www.mrc-bsu.cam.ac.uk/bugs. Accessed Aug. 5, 2006.
Strande´n, I., and D. Gianola. 1999. Mixed effects linear models with
t-distributions for quantitative genetic analysis: A Bayesian approach. Genet. Sel. Evol. 31:25–42.
Wright, S. 1934. An analysis of variability in number of digits in an
inbred strain of guinea pigs. Genetics 19:506–536.
Zellner, A. 1976. Bayesian and non-Bayesian analysis of regression
model with multivariate Student’s t error terms. J. Am. Stat.
Assoc. 71:400–405.
Downloaded from www.journalofanimalscience.org by guest on September 19, 2014
References
This article cites 15 articles, 2 of which you can access for free at:
http://www.journalofanimalscience.org/content/85/6/1363#BIBL
Citations
This article has been cited by 3 HighWire-hosted articles:
http://www.journalofanimalscience.org/content/85/6/1363#otherarticles
Downloaded from www.journalofanimalscience.org by guest on September 19, 2014