Bayesian inference in threshold models using Gibbs sampling

Inference bayesienne dans les modeles a seuil avec echantillonnage de Gibbs. Une analyse bayesienne du modele a seuil avec des categories multiples ordonnees est presentee ici. Les marginalisations necessaires sont obtenues par echantillonnage de Gibbs. On montre que l'utilisation de donnees augmentees - la variable continue sousjacente non observee etant alors consideree comme une inconnue dans le modele conduit a des distributions conditionnelles a posteriori faciles a echantillonner. Celles-ci sont des distributions uniformes independantes pour les seuils et des distributions normales tronquees independantes pour les sensibilites (les variables sous-jacentes). Les parametres restants du modele ont des distributions conditionnelles a posteriori identiques a celles qu'on trouve en modele lineaire gaussien. La methodologie est illustree sur un modele paternel appliquee a une dysplasie de la hanche chez le chien, et les resultats sont compares a ceux d'une etude precedente basee sur un maximum de vraisemblance approche. Deux sequences de Gibbs independantes, longues chacune de 620 000 echantillons, ont ete realisees. Les erreurs d'echantillonnage de type Monte Carlo des moments des densites a posteriori ont ete obtenues par des methodes de series temporelles. Les resultats obtenus avec les 2 sequences independantes sont dans la limite des erreurs d'echantillonnage de Monte-Carlo. A l'exception de la variance paternelle et de l'heritabilite, les distributions marginales a posteriori semblent normales. De ce fait, les inferences basees sur la presente methode sont en bon accord avec celles du maximum de vraisemblance approche. Pour l'estimation des seuils, les sequences de Gibbs revelent de fortes autocorrelations, auxquelles il est cependant possible de remedier en utilisant un autre parametrage.


INTRODUCTION
Many traits in animal and plant breeding that are postulated to be continuously inherited are categorically scored, such as survival and conformation scores, degree of calving difficulty, number of piglets born dead and resistance to disease. An appealing model for genetic analysis of categorical data is based on the threshold liability concept, first used by Wright (1934) in studies of the number of digits in guinea pigs, and by Bliss (1935) in toxicology experiments. In the threshold model, it is postulated that there exists a latent or underlying variable (liability) which has a continuous distribution. A response in a given category is observed, if the actual value of liability falls between the thresholds defining the appropriate category. The probability distribution of responses in a given population depends on the position of its mean liability with respect to the fixed thresholds. Applications of this model in animal breeding can be found in Robertson and Lerner (1949), Dempster and Lerner (1950) and Gianola (1982), and in Falconer (1965), Morton and McLean (1974) and Curnow and Smith (1975), in human genetics and susceptibility to disease. Important issues in quantitative genetics and animal breeding include drawing inferences about (i) genetic and environmental variances and covariances in populations; (ii) liability values of groups of individuals and candidates for genetic selection; and (iii) prediction and evaluation of response to selection. Gianola and Foulley (1983) used Bayesian methods to derive estimating equations for (ii) above, assuming known variances. Harville and Mee (1984) proposed an approximate method for variance component estimation, and generalizations to several polygenic binary traits having a joint distribution were presented by Foulley et al (1987). In these methods inferences about dispersion parameters were based on the mode of their joint posterior distribution, after integration of location parameters. This involved the use of a normal approximation which, seemingly, does not behave well in sparse contingency tables (H6schele et al, 1987). These authors found that estimates of genetic parameters were biased when the number of observations per combination of fixed and random levels in the model was smaller than 2, and suggested that this may be caused by inadequacy of the normal approximation. This problem can render the method less useful for situations where the number of rows in a contingency table is equal to the number of individuals. A data structure such as this often arises in animal breeding, and is referred to as the 'animal model' (Quaas and Pollak, 1980). Anderson and Aitkin (1985) proposed a maximum likelihood estimator of variance component for a binary threshold model.
In order to construct the likelihood, integration of the random effects was achieved using univariate Gaussian quadrature. This procedure cannot be used when the random effects are correlated, such as in genetics. Here, multiple integrals of high dimension would need to be calculated, which is unfeasible even in data sets with only 50 genetically related individuals. In animal breeding, a data set may contain thousands of individuals that are correlated to different degrees, and some of these may be inbred.
Recent reviews of statistical issues arising in the analysis of discrete data in animal breeding can be found in Foulley et al (1990) and Foulley and Manfredi (1991). Foulley (1993) gave approximate formulae for one-generation predictions of response to selection by truncation for binary traits based on a simple threshold model. However, there are no methods described in the literature for drawing inferences about genetic change due to selection for categorical traits in the context of threshold models. Phenotypic trends due to selection can be reported in terms of changes in the frequency of affected individuals. Unfortunately, due to the nonlinear relationship between phenotype and genotype, phenotypic changes do not translate directly into additive genetic changes, or, in other words, to response to selection. Here we point out that inferences about realized selection response for categorical traits can be drawn by extending results for the linear model described .
With the advent of Monte-Carlo methods for numerical integration such as Gibbs sampling (Geman and Geman, 1984;Gelfand et al, 1990), analytical approximations to posterior distributions can be avoided, and a simulation-based approach to Bayesian inference about quantitative genetic parameters is now possible. In animal breeding, Bayesian methods using the Gibbs sampler were applied in Gaussian models by Wang et al (1993Wang et al ( , 1994a   for (co)variance component estimation and by Sorensen et al (1994) and Wang et al (1994b) for assessing response to selection. Recently, a Gibbs sampler was implemented for binary data (Zeger and Karim, 1991) and an analysis of multiple threshold models was described by Albert and Chib (1993). Zeger and Karim (1991) constructed the Gibbs sampler using rejection sampling techniques (Ripley, 1987), while Albert and Chib (1993) used it in conjunction with data augmentation, which leads to a computationally simpler strategy. The purpose of this paper is to describe a Gibbs sample for inferences in threshold models in a quantitative genetic context. First, the Bayesian threshold model is presented, and all conditional posterior distributions needed for running the Gibbs sampler are given in closed form.
Secondly, a quantitative genetic analysis of hip dysplasia in German shepherds is presented as an illustration, and 2 different parameterizations of the model leading to alternative Gibbs sampling schemes are described.

MODEL FOR BINARY RESPONSES
At the phenotypic level, a Bernoulli random variable Y i is observed for each individual i (i = 1, 2, ... , n) taking values y i = 1 or y 2 = 0 (eg, alive or dead).
The variable Y is the expression of an underlying continuous random variable U i , the liability of individual i. When U Z exceeds an unknown fixed threshold t, then Y = 1, and Y = 0 otherwise. We assume that liability is normally distributed, with the mean value indexed by a parameter 0, and, without loss of generality, that it has unit variance (Curnow and Smith, 1975). Hence: where 0' = (b', a') is a vector of parameters with p fixed effects (b) and q random additive genetic values (a), and w' is a row incidence vector linking e to the ith observation.
It is important to note that conditionally on 0, the U i are independent, so for the vector U = {U i } given 0, we have as joint density: where !U(.) is a normal density with parameters as indicated in the argument. In !2!, put WO = Xb + Za, where X and Z are known incidence matrices of order n by p and n by q, respectively, and, without loss of generality, X is assumed to have full column rank. Given the model, we have: where <p(.) is the cumulative distribution function of a standardized normal variate.
Without loss of generality, and provided that there is a constant term in the model, t can be set to 0, and [3] reduces to Conditionally on both 0 and on Y i = y 2 , U i follows a truncated normal distribution. That is, for Yi = 1: where I(X E A) is the indicator function that takes the value 1 if the random variable X is contained in the set A, and 0 otherwise. For Y i = 0, the density is § u ; (w § e, I) /V(-w § e) I (U i x 0) . Invoking an infinitesimal model (Bulmer, 1971), the conditional distribution of additive genetic values given the additive genetic variance in the conceptual base population ( J fl) is multivariate normal: where A is a q by q matrix of additive genetic relationships. Note that a can include animals without phenotypic scores.
We discuss next the Bayesian inputs of the model. The vector of fixed effects b will be assumed to follow a priori the improper uniform distribution: For a description of uncertainty about the additive genetic variance, or a 2, an inverted gamma distribution can be invoked, with density: where v and S' are parameters. When v = -2 and S Z = 0, [8] reduces to the improper uniform prior distribution. A proper uniform prior distribution for Q a is: where k a is a constant and a a 2m!'. is the maximum value which J£ can take a priori.
To facilitate the development of the Gibbs sampler, the unobserved liability U is included as an unknown parameter in the model. This approach, known as data augmentation (Tanner and Wong, 1987;Gelfand et al, 1992;Albert and Chib, 1993;Smith and Roberts, 1993) leads to identifiable conditional posterior distributions, as shown in the next section.
Bayes theorem gives as joint posterior distribution of the parameters: The last term is the conditional distribution of the data given the parameters. We notice that, for Y i = 1, say, we have For Y i = 0, we have: This distribution is degenerate, as noted by Gelfand et al (1992) because knowledge of U i implies exact knowlege of Y i . This can be written (eg, Albert and Chib, 1993) as: The joint posterior distribution [10] can then be written as: where the conditioning on hyperparameters v and S' is replaced by 0 a ma . when the uniform prior [9] for the additive genetic variance is employed.

Conditional posterior distributions
In order to implement the Gibbs sampler, all conditional posterior distributions of the parameters of the model are needed. The starting point is the full posterior distribution !13!. Among the 4 terms in (13!, the third is the only one that is a function of b and we therefore have for the fixed effects: which is proportional to !U(Xb+Za, I). As shown in Wang et al (1994a), the scalar form of the Gibbs sampler for the ith fixed effect consists of sampling from: where x i is the ith column of the matrix X, and b i satisfies: In !16!, X-i is the matrix X with the column associated with i deleted, and b_ i is b with the ith element deleted. The conditional posterior distribution of the vector of breeding values is proportional to the product of the second and third terms in !13!: which has the form !(0,Acr!)!(u!b,a). Wang et al (1994a) showed that the scalar Gibbs sampler draws samples from: where z i is the ith column of Z, c ii is the element in the ith row and column of A-1 , /B B = (Qa)-1, and a i satisfies: In [19], c i ,-i is the row of A-1 corresponding to the ith individual with the ith element excluded. We notice from [14] and [17], that augmenting with the underlying variable U, leads to an implementation of the Gibbs sampler which is the same as for the linear model, with the underlying variable replacing the observed data.
For the variance component, we have from !13!: Assuming that the prior for o,2is the inverted gamma given in !8!, this becomes: and assuming the uniform prior !9!, it becomes: Expression !21a! is in the form of a scaled inverted gamma density, and [21b] in the form of a truncated scaled inverted gamma density.
The conditional posterior distribution of the underlying variable U i is proportional to the last 2 terms in !13!. This can be seen to be a truncated normal distribution, on the left if Y i = 1 and on the right otherwise. The density function of this truncated normal distribution is given in !5!. Thus, depending on the observed Yi, we have: or Sampling from the truncated distribution can be done by generating from the untruncated distribution and retaining those values which fall in the constraint region. Alternatively and more efficiently, suppose that U is truncated and defined in the interval !i, j] only, where i and j are the lower and upper bounds, respectively. Let the distribution function of U be F, and let v be a uniform [0, 1] variate. Then U = F-1 !F(i) + v(F(j) &mdash; F(i))! is a drawing from the truncated random variable (Devroye, 1986). Albert and Chib (1993) also constructed the mixed model in terms of a hierarchical model, but proposed a block sampling strategy for the parameters in the underlying scale, instead. Essentially, they suggest sampling from the distributions ( Q a, a, blU) as (a!IU)(a, blU, J£ ), instead of from the full conditional posterior distributions !15!, !18! and !21!, and they assumed a uniform prior for log(a2) in a finite interval. To facilitate sampling from p(or2lU), they use an approximation which consists of placing all prior probabilities on a grid of or2 values, thus making the prior and the posterior discrete. The need for this approximation is questionable, since the full conditional posterior distribution of (T has a simple form as noted in [21] above. In addition, in animal breeding, the distribution (a, b[U, a a 2) is a high dimensional multivariate normal and it would not be simple computationally to draw a large number of samples.

MULTIPLE ORDERED CATEGORIES
Suppose now that the observed random variable Y can take values in one of C mutually exclusive ordered categories delimited by C + 1 thresholds. Let to = -oo, t c = +oo, with the remaining thresholds satisfying t1 ! t 2 ... < t C -1 -

Generalizing [3]:
Conditionally on A, Y i = j, t j -1 and t j , the underlying variable associated with the ith observation follows a truncated normal distribution with density: Assuming that o, a, 2 b and t are independently distributed a priori, the joint posterior density is written as: where p(Ulb, a, t) = p(Ulb, a). Generalizing [12], the last term in [25] can be expressed as (Albert and Chib, 1993): All the conditional posterior distributions needed to implement the Gibbs sampler can be derived from !25!. It is clear that the conditional posterior distributions of b i , a i and u2 are the same as for the binary response model and given in (15!, [18] and !21!. For the underlying variable associated with the ith observation we have from !25!: This is a truncated normal, with density function as in !24!. The thresholds t = (t l , t 2 , ... , tC-1) are clearly dependent a priori, since the model postulates that these are distributed as order statistics from a uniform distribution in the interval [t mini t max] . However, the full conditional posterior distributions of the thresholds are independent. That is, p(t j It-j , b, a, U, o l a 2 , Y ) = p(t j I U, y), as the following argument shows. The joint prior density of t is: where T = {(h, t2,&dquo;', tc-dltmin ! / x t2 ! ... ! t C -1 v t max ) (Mood et a l , 1974). Note that the thresholds enter only in defining the support of p(t). The conditional posterior distribution of t j is given by: which has the same form as !26!. Regarded as a function of t, [26] shows that, given U and y, the upper bound of threshold t j is min (U I Y = j +1) and the lower bound is max(UIY = j). The a priori condition t E T is automatically fulfilled, and the bounds are unaffected by knowledge of the remaining thresholds. Thus t j has a uniform distribution in this interval given by: This argument assumes that there are no categories with missing observations. To accommodate for the possibility of missing observations in 1 or more categories, Albert and Chib (1993) define the upper and lower bounds of threshold j, as minfmin(UIY = j + 1), t j+d and as max{max(U!Y = j),t j -' 11 , respectively. In this case, the thresholds are not conditionally independent. The Gibbs sampler is implemented by sampling repeatedly from !15!, !18!, !21!, [24] and (28!. Alternative parameterization of the multiple threshold model The multiple threshold model can also be parameterized such that the conditional distribution of the underlying variable U, given 0, has unknown variance 0' instead of unit variance. The equivalence of the 2 parameterizations is shown in the Appendix. This parameterization requires that records fall in at least 3 mutually exclusive ordered categories; for C categories, only C-3 thresholds are identifiable.  (1988). Briefly, the records consisted of radiographs of 2 674 offspring from 82 sires. These radiographs had been classified according to guidelines approved by FCI (Federation Cynologique Internationale, 1983), each offspring record was allocated to 1 of 7 mutually exclusive ordered categories. The model for the underlying variable was: . where a i is the effect of sire i (i = 1, 2, ... , 82; j = 1, 2,... , n i ). The prior distribution of /! was as in [7] and sire effects were assumed to follow the normal distribution: The prior distribution of the sire variance ( Q a) was in the form given in !8!, with v = 1 and S 2 = 0.05. The prior for t 2 , ... , t6 was chosen to be uniform on the ordered subset of [f, = -1.365, f 7 = +00!5 for which ti < t 2 < ... < t 7 , where f l was the value at which t 1 was set, and f 7 is the value of the 7th threshold.
The value for f l was obtained from Andersen et al (1988), in order to facilitate comparisons with the present analysis. The analysis was also carried out under the parameterization where the conditional distribution of U given 0 has variance a 2 Here, Q e was assumed to follow a prior of the form of !8!, with v = 1 and S Z = 0.05 and t 6 was set to 0.429. Results of the 2 analyses were similar, so only those from the second parameterization are presented here.

Gibbs sampler and post Gibbs analysis
The Gibbs sampler was run as a single chain. Two independent chains of length 620 000 each were run, and in both cases, the first 20 000 samples were discarded. Thereafter, samples were saved every 20 iterations, so that the total number of samples kept was 30 000 from each chain. Start values for the parameters were, for the case of chain 1, o,2= 2.0, Q a = 0.5, t 2 = -0.8, t 3 = -0.5, t 4 = -0.2, t 5 = 0.1. For chain 2, estimates from Andersen et al (1988) were used, and these were J fl = 1.0, or = 0.1, t 2 = -1.05, t 3 = -0.92, t 4 = -0.62, t 5 = -0.34. In both runs, starting values for sire effects were set to zero. Two important issues are the assessment of convergence of the Gibbs sampler, and the Monte-Carlo error of estimates of features of posterior distributions. Both issues are related to the question of whether the chain, or chains, have been run long enough. This is an area of active research in which some guidelines based on theoretical work (Roberts, 1992;Besag and Green, 1993;Smith and Roberts, 1993;Roberts and Polson, 1994) and on practical considerations (Gelfand et al, 1990;Gelman and Rubin, 1992;Geweke, 1992;Raftery and Lewis, 1992) have been suggested. The approach chosen here is based on Geyer (1992), who used time series methods to estimate the Monte-Carlo error of moments estimated from the Gibbs chain. Other approaches include, for example, batching (Ripley, 1987), and Raftery and Lewis (1992) proposed a method based on 2-state Markov chains to calculate the number of iterations needed to estimate posterior quantiles. Let X l , ... , X&dquo;, be elements of a single Gibbs chain of length m. The m sampled values, which are generally correlated, can be used to compute features of the posterior distribution. In our model, the Xs could represent samples from the posterior distribution of a sire value, or of the sire variance, or of functions of these. Invoking the ergodic theorem, for large m, the posterior mean can be estimated by and the variance of the posterior distribution can be estimated by (Cox and Miller, 1965;Geyer, 1992). Similarly, a marginal posterior distribution function, F(Z) = P(X < Z), can be estimated by which means that F(Z) is estimated by the empirical distribution function. All these estimators are subject to Monte-Carlo sampling error, which is reduced by prolongation of the chain.
Consider the sequence g(X 1 ),..., g(Xm), where g(.) is a suitable function, eg, g(X) = X or g(X) = l!X!Z!, and let E(g(X)) = p. The lag-time auto-covariance of the sequence is estimated as: where is the sample mean for a chain of length m, and t is the lag. The auto-correlation is then estimated as Tm(!)/!m(0). The variance of the sample mean of the chain is given by which will exceed y(0)/m if q(t) > 0 for all t, as is usual. Several estimators of the variance of the sample mean have been proposed (Priestley, 1981), but we chose one suggested by Geyer (1992), which he calls the initial positive sequence estimator. Let F m (t) = $ m(2t) + f im (2t + 1), t = 0, 1, .... The estimator can then be written as where t is chosen such that it is the largest integer satisfying fm (I) > 0, i = 1, 1 ... t.
The justification for this choice is that r(i) is a strictly positive, strictly'decreasing function of i. If X l , X 2 , ... X&dquo;, are independent, then Var(í 1m ) = q(0) /m. To obtain an indication of the effect of the correlation on Var(í 1m ), an 'effective number' of independent observations can be assessed as j m = 5 m (0)/ Vâr(í 1m ) ' When the elements of the Gibbs chain are independent, !m = m.

RESULTS
Estimates of the empirical distribution function of various parameters of the model for each of the 2 chains of the Gibbs sampler are shown in figures 1-5. For example, figure 2 shows that there is a 90% posterior probability that the sire variance lies between 0.065 and 0.14, and the median of this posterior distribution is slightly under 0.10. Similarly, figure 3 indicates that there is 90% posterior probability that heritability in the underlying scale (h 2 = 4J£ / (J£ + Jfl ) ) lies between 0.24 and 0.49, and the median of the posterior distribution is 0.35. Although this distribution is slightly skewed, the estimate of the median agrees well with the ML type estimate of heritability of 0.35, reported in Andersen et al (1988). Figure 4 depicts estimates of distribution functions for the mean ( J1 ) and for each of 3 sire effects (a l , a 2 , a 3 ). Figure 5 gives corresponding distributions for 2 threshold parameters (t 2 , and t 3 ).
The figures fall in 3 categories. The distribution functions obtained from chains 1 and 2 coincide for each of the variables Q a, h2, a l , a 2 and a 3 , where the sire effects a l , a 2 and a 3 pertain to 3 males with 31, 5 and 158 offspring, respectively. A small deviation between chains 1 and 2 is observed for Jfl and J1 , and a larger deviation is observed for the threshold parameters ( fig 5). The Gibbs sequence for the threshold parameters showed very slow mixing properties. For example, for threshold 2, the autocorrelations between sampled values were 0.785, 0.663 and 0.315, for lags between 5, 10 and 50 samples, respectively. The reason is that the sampled value for a given threshold is bounded by the values of the neighbouring underlying variables U. If these are very close, the value of the threshold in subsequent samples is likely to change very slowly. Under the parameterization where 1 of the thresholds is substituted by the residual variance, the autocorrelations associated with the lags above, between samples from the marginal posterior distribution of e 2, were 0.078, 0.064 and 0.032, respectively. Another scheme that may accelerate mixing is to sample jointly from the threshold and the liability. For sire effects, lag autocorrelations were close to zero.
A comparison between marginal posterior means (average of the 30 000 samples) estimated from the 2 chains is shown in table I. The difference between chains 1 and 2 is in all cases within Monte-Carlo sampling error, which was estimated within chains, using [34]. The 'effective number' of observations u Jm for the means of marginal posterior distributions is close to 30 000 for or a 2, a l , a 2 and a 3 , but is about 3 000 for Q e and p and between 200 and 400 for t 2 , ... , ts.
The marginal posterior distributions for J1 , t 2 , ... , t 5 , a l , a 2 , a 3 are well approximated by normal distributions. The posterior means and standard deviations of these marginal distributions can be compared to estimates reported in Andersen et al (1988), who used a 2-stage procedure. The authors first estimated variances using an REML-type estimator (Harville and Mee, 1984). Secondly, assuming that the estimated variances were the true ones, fixed effects and sire effects were estimated as suggested by Gianola and Foulley (1983). For example, for the 3 sires with 158, 31 and 5 offspring respectively, the 2-stage procedure yielded estimates of sire effects (approximate posterior standard deviations) of 0.30 (0.09), 0.17 (0.17) and -0.093 (0.26), respectively. The present Bayesian approach with the Gibbs sampler, yielded estimates of marginal posterior means and standard deviations for these sires of 0.304 (0.088), 0.177 (0.166), and -0.092 (0.263), respectively (table I).

DISCUSSION
We have described a Gibbs sampler for making inferences from threshold models for discrete data. The method was illustrated using a model with unrelated sires; here likelihood inference with numerical integration is a competing alternative. For this model and data, marginal posterior distributions of sire effects are well approximated by normal distributions. On the other hand, with little information per random effect, eg, animal models, the normality of marginal posterior distributions when variances are not known is unlikely to hold. A strength of the Bayesian approach via Gibbs sampling is that inferences can be made from marginal posterior distributions in small sample problems, without resorting to asymptotic approximations. Further, the Gibbs sampler can accommodate multivariately distributed random effects, such as is the case with animal models, and this cannot be implemented with numerical integration techniques. It seems important to investigate threshold models further, especially with sparse data structures consisting of many fixed effects with few observations per subclass. This case was studied by Moreno et al (manuscript submitted for publication) in the binary model, where they investigated frequency properties of Bayesian point estimators in a simulation study. They showed that improper uniform prior distributions for the fixed effects lead to an average of the marginal posterior mean of heritability which was larger than the simulated value. They obtained better agreement when fixed effects were assigned proper normal prior distributions. The use of 'non-informative' improper prior distributions is discouraged on several grounds by, among others, Berger and Bernardo (1992), as this can lead to improper posterior distributions. It seems that the disagreement between simulated and estimated values in Moreno et al is due to lack of information and not due to impropriety of posterior distributions. Thus, the bias persists, though smaller, when all parameters of the model are assigned proper prior distributions (Moreno, personal communication). It is clear that as the number of fixed effects increases, for a constant amount of observations, a larger proportion of fixed effect levels will contain data falling into only one of the 2 dichotomies. There is no information in the data (in the Fisherian sense) to estimate these fixed effects and the likelihood is ill-conditioned. In the case of these sparse data structures, the choice of the prior distribution for the fixed effects may well be the most critical part of the problem. This needs to be studied further.
Data augmentation in the Gibbs sampler led to conditional posterior distributions which are easy to sample from. This facilitates programming. We have noted though, that threshold parameters have very slow mixing properties, and this is probably related to the data augmentation approach used in this study (Liu et al, 1994). With our data, the parameterization in terms of the residual variance resulted in smaller autocorrelations between samples of Q e than between samples of the thresholds. A scheme that is likely to accelerate mixing is to sample jointly from the threshold and liability. This step may necessitate other Monte-Carlo sampling techniques such as a Metropolis algorithm (Tanner, 1993), since sampling is from a non-standard distribution. Alternative computational strategies and parameterizations of the model may be more critical with animal models. Here, there is typically little information on additive genetic effects, and these are correlated.
These properties slow down convergence of the Gibbs chain.
The methods described in this paper can be adapted easily to draw inferences about genetic change when selection is for categorical data. Sorensen et al (1994) described how to make inferences about response to selection in the context of the Gaussian model. In the threshold model, the only difference is that observed data are replaced by the unobserved underlying variable U (liability). In order to make inferences about response to selection, the parameterization must be in terms of an animal model.