On biased inferences about variance components in the binary threshold model

Une etude de simulation a ete conduite pour etudier les proprietes frequentistes de trois estimateurs de la composante de variance dans un modele mixte a seuils pour variable binaire. Les trois estimateurs ont ete : le mode de l'approximation normale de la distribution marginale a posteriori de la composante qui est appelee dans la litterature « maximum de vraisemblance marginale» (MML), l'esperance de distribution marginale a posteriori de la composante, en utilisant l'echantillonnage de Gibbs (GSR) et finalement, le mode de la distribution conjointe a posteriori des parametres de position et de variance, utilise avec une methode iterative de correction des biais par bootstrap (MJP-IBC). Cette derniere a ete recemment proposee dans la litterature comme une methode d'obtention d'estimateurs approximativement non biaises. Les resultats de l'etude confirment que le MML peut donner des inferences biaisees a propos de la composante de variance et que le signe du biais depend de la quantite d'information attachee aux effets fixes ou aux effets aleatoires. Le GSR peut produire des estimees a biais positif quand il y a peu d'information par niveau d'effet fixe. Dans ce cas, le biais persiste meme quand l'information sur la composante de variance est importante et quand la distribution posterieure est propre. La distribution marginale a posteriori est hautement concentree et symetrique mais montre un decalage a droite de la vraie valeur (simulee). Ce biais peut etre reduit en considerant que la distribution a priori des effets fixes est de type gaussien mais cette strategie n'est pas efficace pour les structures de donnees ou l'information est tres dispersee. La methode MJB-IBC a produit des estimees non biaisees de la composante de variance dans tous les cas etudies. Cet estimateur est simple a calculer mais contrairement a GSR avec des effets fixes a distribution gaussienne, elle peut mener a des estimees situees en dehors de l'espace des parametres.

Summary -A simulation study was conducted to study frequentist properties of three estimators of the variance component in a mixed effect binary threshold model. The three estimators were: the mode of a normal approximation to the marginal posterior distribution of the component, which is denoted in the literature as marginal maximum likelihood (MML); the mean of the marginal posterior distribution of the component, using the Gibbs sampler to perform the marginalisations (GSR); and third, the mode of the joint posterior distributions of location and the variance parameter, used in conjunction with the iterative bootstrap bias correction (MJP-IBC). The latter was recently proposed in the literature as a method to obtain nearly unbiased estimators. The results of this study confirm that MML can yield biased inferences about the variance component, and that the sign of the bias depends on the amount of information associated with either fixed effects or with random effects. GSR can produce positively biased inferences when the amount of data per fixed effect is small. When fixed effects are poorly estimated, the bias persists, despite the fact that posterior distributions are guaranteed to be proper, and that the amount of information about the variance component is large. In this case, the marginal posterior distribution of the variance component is highly peaked and symmetric, but it shows a shift towards the right with respect to the true (simulated) value. This bias can be reduced by assigning a Gaussian probability density function to the prior distribution of the fixed effects, but this strategy does not work with very sparse data structures. The method based on MJP-IBC yielded unbiased inferences about the variance component in all the cases studied. This estimator is computationally simple, but contrary to GSR with normal priors for the fixed effects, can lead to estimates that fall outside the parameter space.
binary traits / Gibbs sampling / biased estimators / bootstrap * Correspondence and reprints Résumé -À propos des biais d'estimation de composantes de variance dans le modèle binaire à seuil. Une étude de simulation a été conduite pour étudier les propriétés fréquentistes de trois estimateurs de la composante de variance dans un modèle mixte à seuils pour variable binaire. Les trois estimateurs ont été : le mode de l'approximation normale de la distribution marginale a posteriori de la composante qui est appelée dans la littérature « maximum de vraisemblance marginale » (MML), l'espérance de distribution marginale a posteriori de la composante, en utilisant l'échantillonnage de Gibbs (GSR) et finalement, le mode de la distribution conjointe a posteriori des paramètres de position et de variance, utilisé avec une méthode itérative de correction des biais par bootstrap (MJP-IBC). Cette dernière a été récemment proposée dans la littérature comme une méthode d'obtention d'estimateurs approximativement non biaisés. Les résultats de l'étude confirment que le MML peut donner des inférences biaisées à propos de la composante de variance et que le signe du biais dépend de la quantité d'information attachée aux effets fixes ou aux effets aléatoires. Le GSR peut produire des estimées à biais positif quand il y a peu d'information par niveau d'effet fixe. Dans ce cas, le biais persiste même quand l'information sur la composante de variance est importante et quand la distribution postérieure est propre. La distribution marginale a posteriori est hautement concentrée et symétrique mais montre un décalage à droite de la vraie valeur (simulée). Ce biais peut être réduit en considérant que la distribution a priori des effets fixes est de type gaussien mais cette stratégie n'est pas e,!cace pour les structures de données où l'information est très dispersée. La méthode MJB-IBC a produit des estimées non biaisées de la composante de variance dans tous les cas étudiés. Cet estimateur est simple à calculer mais contrairement INTRODUCTION Ordered categorical traits are often analysed by means of the threshold liability model, 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. Here we will consider the situation where the liability is modelled as a linear function of fixed effects and possibly correlated random effects. In a typical likelihood setting, one may be interested in making joint inferences about estimable functions of fixed effects and on variance components associated with the random effects. Within a Bayesian framework, one may concentrate inferences on the marginal posterior distributions of contrasts of 'fixed' effects and of the variance components, and possibly on the marginal posterior distributions of the random effects. In a genetic context, the latter involve typically additive genetic values and knowledge about these is important to carry out selection decisions and to describe response to selection. Other inferences of interest may require computing probabilities that observations fall in a given category for given individuals or combinations of parameters.
All these inferential problems, whether in a likelihood or in a Bayesian setting, require marginalisations. For example, computation of the likelihood or of the joint posterior distribution of fixed effects and variance components requires integration with respect to the random effects. Although these integrals do not have a closed form solution, numerical integration using quadrature is feasible when the random effects are uncorrelated (Anderson and Aitkin, 1985). When random effects are correlated, as is the case with animal models, these numerical integration techniques cannot be used. To circumvent this problem, a large number of approximations and different approaches has been suggested in the literature (see, for example, McCullagh and Nelder (1989) for a discussion of some of these, and Foulley et al (1990) and Foulley and Manfredi (1991) for a review of the analysis of categorical traits in an animal breeding context).
The introduction of Markov chain Monte Carlo methods allows computation of multidimensional integrals and analytic approximations can therefore be avoided. The use of the Gibbs sampler to analyse ordered categorical traits has been reported by Zeger and Karim (1991) and by Albert and Chib (1993), using a Bayesian perspective. McCulloch (1994) applies the Gibbs sampler to estimate variance components for binary data in a likelihood setting. In an animal breeding context, the Gibbs sampler was applied to the analysis of categorical traits by Moreno (1993) and by Sorensen et al (1995), both from a Bayesian perspective.
Zeger et al (1988) studied non-linear models for the analysis of longitudinal data, and drew attention to the dependence between the marginal expectation of the data and the variance of the random effects. In the context of a generalized estimating equations approach, they showed that estimators of the fixed effects and of the variance component associated with the random effects are not even asymptotically orthogonal, as they are in the Gaussian model. A similar dependency in the case of the logistic mixed model was shown by Drum and McCullagh (1993) in the context of restricted maximum likelihood, where the latter was constructed using error contrasts. This approach to constructing the restricted likelihood generates an equation that is not a function of the fixed effects in the case of the Gaussian model, but not in the case of non-linear models. The association between the fixed effects and the variance component indicates that inference about the latter is dependent upon the amount of information in the data with which fixed effects are estimated.
The study of Moreno (1993) reported frequentist properties of the mean of the marginal posterior distribution of the variance component in a binary mixed threshold model. In this study, a sire model was assumed. Moreno (1993) noted that when there is little information about fixed effects, the mean of the marginal posterior distribution of the variance component is biased upwards. This is so even in cases where there are very large numbers of sires and offspring per sire, and consequently this marginal posterior distribution is peaked and symmetric. It is well known that when the data associated with a particular fixed effect fall all within the same category (known as the extreme category problem (ECP), the likelihood is underidentified (eg, Misztal et al, 1989;Gelman et al, 1995). In this situation, in a Bayesian setting, the assignment of improper prior distributions to the 'fixed effects' can lead to improper posterior distributions. However, we show that the bias related to inference about the variance component when there is little information about fixed effects persists, when all the parameters of the model are assigned proper priors. In this case, the posterior distribution is guaranteed to be proper.
This bias problem has been the subject of an extensive simulation study by Hoeschele and Tier (1995). These authors, who used the Gibbs sampler with data augmentation (Tanner and Wong, 1987), confirm the results of Moreno (1993), and indicate that part of the problem can be alleviated by assuming that 'fixed effects' follow a priori a Gaussian distribution. As we show below, treating 'fixed effects' as if they were 'random' can reduce the bias associated with the mean of the marginal posterior distribution of the variance component, but this is not a strategy that always works. An alternative approach to inferences in the mixed threshold model was presented by Tempelman (1994). The method is quite general and is based on an approximation to the marginal posterior distribution of the variance component, where the marginalisation is obtained using Laplacian integration. The computational burden of this method amounts to obtaining repeatedly the determinant of the Hessian of the joint posterior distribution, conditional on the current value of the variance component. The dimension of the Hessian is equal to the number of 'fixed' and random effects in the model. In a small simulation study, Tempelman (1994) showed that the method performed well in a sire model, but yielded biased estimates of the variance component in an animal model. Breslow and Lin (1995) have recently developed methods to adjust for the bias in approximate estimators of the variance component and regression coefficients for generalised linear mixed models having one source of random variation. The bias correction applies to approximations based on penalised quasilikelihood and holds for generalised linear mixed models having canonical link functions, the latter being a condition that is not satisfied in the case of the probit mixed model studied in this work. Kuk (1995) proposed an alternative iterative method using the bootstrap to obtain asymptotically unbiased estimators. The appealing aspect of this approach is that in principle it can be applied to any estimator. Though computationally intensive, the method can be implemented in a straightforward manner. The objective of this work is to study frequentist properties of three approaches to making inferences about variance components in mixed binary threshold models. One of these approaches, which does not rely on analytic approximations, is based on the mean of the marginal posterior distribution of the variance components. We study its properties using different data structures and assuming either proper or improper prior distributions for the 'fixed' effects and for the variance component.
We report also frequentist properties of a computationally simple estimator based on the mode of the joint posterior distribution of location and scale parameters, used in conjunction with the iterative bias correction using the bootstrap proposed by Kuk (1995).

MATERIAL AND METHODS
Model for the analysis of binary responses A Bernoulli random variable yg is recorded for observation j(j = 1,...,!,) in subclass i(i = 1, ... , s), and takes the value y2! = 1 if !i! > t or y2! = 0 otherwise, where Q2! represents the underlying variable (liability) and t is the threshold.
Conditionally on fixed effects (herds) (b, of order p) and on random (sire) effects (u, of order q), the !z! are independently and normally distributed, ie, N(x'b + z!u, 1), and the Yij are independent with Pr(yg = O!b, u) _ 4l(t -xib -z!u).
The symbol 4l(.) denotes the standard normal cumulative distribution function and x!(z!) denotes the ith row of the known design matrix X(Z). The subclass i is defined as the herd-sire combination i. We will denote the prior distribution for the fixed effects as p(b). Invoking an infinitesimal model, the vector of sire effects is multivariate normally distributed: ulA, a 2 -N(O, A U2 ), where A is a known q by q matrix af additive genetic relationships among sires, and Q u is the unknown variance between sires. The prior probability density for Q u will be denoted p(o, u 2).
Under the model, the conditional distribution of the data y, given b and u is: where n il is the number of observations in class i with !2! = 0 and n 22 is the number of observations in class i with y2! = 1. The joint posterior distribution is given by: Unless otherwise stated, p(b) is assumed uniform in the interval [t -5, t + 5!, and p( Q u) is assumed uniform in the interval !0,1/3!. Consequently, in these cases, posterior distributions are proper.
As is well established (Cox and Snell, 1989), in the binary threshold model neither the residual variance nor the threshold can be estimated. In this work, the simulated data are analysed using a model in which the residual variance is set equal to 1, and the value of the threshold is arbitrarily set equal to the simulated value.

Methods of inference
In this section we briefly describe the three methods of inference studied in this work. These are, first, marginal maximum likelihood (MML) .
The second method is based on an estimator of the mean of the marginal posterior distribution of the sire variance. Here, marginalisations were obtained using the Gibbs sampler applying an adaptive rejection sampling algorithm due to Gilks and Wild (1992) (GSR). The third method is based on the mode of the joint posterior distribution !2!, used in conjunction with the iterative bootstrap bias correction (MJP-IBC).
The method based on MML was chosen because it is one of the most common methods used in the analysis of categorical traits. Method GSR does not involve analytic approximations and therefore it is interesting to compare it to MML under a variety of data structures. As we show below, both MML and GSR in the form implemented in this work, can lead to biased inferences, even in cases with a large number of sires and many offspring per sire. Method MJP-IBC is an ad hoc method which yields inferences with no detectable bias.

Marginal maximum likelihood (MML)
This method is based on locating the modal value of an approximation to In p(QuIY)! the logarithm of the marginal posterior distribution of the variance component.
As shown  and applying their result for the case Q u u 2 _ Uniform [0, 1/3]: The expression on the right hand side of [3] involves finding the expected value of u'A-l u, where this expectation is taken with respect to the distribution of uly, afl . Var(uly, (7!). In the context of the present binary model, the form of the distribution of uly, (7! is not known. MML consists of approximating E(uly, (7!) by the mode of p(b, uJy, (7!) and Var(uly, (7!) by the inverse of minus the expected value of the second derivatives of p(b, uly, (7!) with respect to b and u. Setting [3] equal to zero and solving for the unknown, leads to an iterative scheme which involves the solution of a non-linear system reminiscent of Henderson's (1973)  where C uu is the part of the complete inverse of the coefficient matrix corresponding to the u by u block.
This method is known to yield biased inference in small samples with respect to the variance component Dellaportas and Smith, 1993;Engel et al, 1995). The sign of the bias depends on the data structure (Hoeschele and Tier, 1995;Engel et al, 1995). With a large amount of information on fixed effects and little information on random effects, the bias is negative. When the amount of information on the random effects is large, but fixed effects are poorly estimated, the bias is positive. This behaviour of the estimator can lead to situations where no bias is detected.
Mean of the marginal posterior distribution of sire variance, using the Gibbs sampler (GSR) The Gibbs sampler is an iterative numerical technique for the approximate generation of (possibly dependent) samples from multivariate distributions. The iterative scheme consists of obtaining draws from the fully conditional posterior distributions of all the parameters of the model. The literature on Markov chain Monte Carlo methods is very large and many theoretical and applied issues can be found in the recent publication of Gilks et al (1996).
In many applications, the fully conditional posterior distributions are known, and the implementation of the Gibbs sampler is straightforward. This is the case of the binary threshold mixed model, when the Gibbs sampler is used in conjuction with data augmentation, as in Albert and Chib (1993). If the technique of data augmentation is not used, the present model leads to fully conditional posterior distributions of the location parameters that are not of standard form. In this case, other sampling schemes can be implemented, such as the adaptive rejection sampling proposed by Gilks and Wild (1992). This is the approach chosen in this paper. Details of the fully conditional posterior distributions can be found in the Appendix. We follow closely the adaptive rejection sampling algorithm as described in Dellaportas and Smith (1993).
The implementation of GSR requires draws from the fully conditional posterior distributions of the parameters. One of these is the fully conditional posterior distribution of the sire variance, p ( 0 , 2 u 1 b, u, y), which is in the form of a scaled inverted chi-square. The heritability, h 2 = 40&dquo;! u / (0,2 + 1) , was constructed for each value drawn from p(a u 2 lb, u, y). The Monte Carlo estimate of the mean of the marginal posterior distribution of heritability was computed as the mean h 2 over the Gibbs samples. In the present study, differences between mean, mode and median of the marginal posterior distribution of heritability were negligible and results based on the mean only are reported. The Gibbs sampler was implemented using several chains, and all the samples in each were kept. Burn in periods were determined on a pragmatic basis. The criterion of Gelman and Rubin (1992) was adopted to determine convergence.
The Gibbs sampler allows the drawing of inferences using a fully Bayesian approach, without resorting to analytical approximations. Asymptotically, posterior distributions are normal, with mean equal to the maximum likelihood estimator.
It seems relevant to study the behaviour of this approach under different data structures, in view of the fact that analytic approximations are not involved.
Mode of the joint posterior distribution followed by the iterative boostrap bias correction (MJP-IBC) The first step of this method is to find the mode of the joint posterior distribution of location and dispersion parameters !2!, under the assumption that the prior of p( 0 ,2) follows a scaled inverted chi-square distribution. The latter is of the form: where q * and u* 2 are parameters of the prior distribution. In the present study, q * was set equal to q, and o,* 2 was derived assuming h 2 = 0.9. Other starting values of h 2 were also tried, but convergence was faster using this value. Differentiating the resulting joint posterior with respect to o, u 2, setting to zero and solving for or2 U ) produces the iterative scheme: as shown by Hoeschele et al (1987), where u is the mode of p(b, ula!,y) as in MML. Notice that in contrast to !5!, [7] does not require knowledge of the inverse of any matrix and is very easy to compute. Use of a proper prior distribution for o l2 is important, because with an improper uniform prior, this method is known to yield zero estimates for Q u (Lindley and Smith, 1972).
The second step, consists of applying the iterative bias correction that was proposed by Kuk (1995). This procedure can be described as follows. Let 0 be the parameter of interest, and let e represent an initial estimate obtained as a solution to a system of estimating equations !(0!y), where y are the data that are assumed to follow a distribution f(y!0). The system of equations !(9!y) could be the score function, for example, or some approximation to it. In the present case, ljJ(ely) is the system of equations that results from the joint maximization of the posterior distribution with respect to location and dispersion parameters. Let A * , the expected value of the estimator of e, satisfy: If the estimator of 0 is biased, the asymptotic bias is given by: where g(e) = e * . Kuk (1995) proposes the following correction for this bias. Let (3!°! be an initial estimate of the bias of 0. Then define an updated estimate of the bias as: ...

I-
1 .1 I-. 11 and define an updated bias-corrected estimate of 0 as: Kuk (1995) shows that this estimator is asymptotically consistent. Since the function g (8) is not usually known, Kuk (1995) proposes to approximate g(8) by g M (8) ! which is the average of e over simulated samples. That is: where y l , ... , y&dquo;, are simulated samples (ie, data) from f (y!9), and i (y i ) is the estimate of e obtained using the simulated sample i, y 2 . Replacing g by g M in !10!, Kuk (1995) obtains: The simulated values y i ,..., y m are obtained from f (y!0), where in the kth round of iteration, e is set to 8 -(3!).
The implementation of the above procedure using method (MJP-IBC) was as follows. First, estimates of location and scale parameters, e, are obtained from the joint maximization of the posterior distribution !2!, conditional on the original data. Second, simulated samples are generated from f(ylO), where e is set to the value of these initial estimates 6. Third, 6(y i ) are obtained from the joint maximization of the posterior distribution [2], conditional on the ith simulated sample (i = 1, ... , m). Fourth, an average of 8( Yi ) over the m samples is computed, which is denoted by g M (9), and the estimate of the bias is obtained from [13]. Finally, 6 is corrected using !14!, which is used to generate a new set of simulated samples. This iterative scheme is repeated until convergence is achieved.
In the presence of the extreme category problem (ECP), it was required that the proportion of ECPs I n the simulated samples was the same as in the original data. This was achieved by adjusting the values of the fixed effects associated with ECPs. Another point to highlight in the application of MJP-IBC in models with ECPs is that the correction in expression [14] involves the component of variance and the non-ECP fixed effects only. Moreover, if in a simulated sample y i , an ECP is generated for a fixed effect that was a non-ECP fixed effect in the original data, then the element in 6(y i ) associated with this fixed effect, is set to the value that ' was used to simulate this fixed effect in the sample y i , and not to its estimated value. This ad hoc strategy was arrived at on a trial and error basis, and it was followed because it avoided problems of lack of convergence. Other strategies could be envisaged, one such being based on assigning mildly informative priors to the 'fixed' effects.

Simulation study
A Monte Carlo simulation study was carried out in order to evaluate the three methods of inference under a variety of data structures. The criteria used to assess the performance of the methods were bias and mean square error. The data were simulated as follows. Liability was drawn from a normal distribution, with mean equal to the fixed effect corresponding to the individual in question and variance given by (72 + o,,2, where U2 is the sire variance and U2 was set equal to 1. Only one set of fixed effects were simulated (herds), and these were drawn from uniform distributions within the interval (-2, 2). Sires were drawn from a normal distribution, with mean (1/2)v,P and variance (3/4) Q u, where up is the realized value of the sire of the individual. Base population sires were drawn assuming a mean of zero and variance u U. 2 Three non-overlapping generations of data (with no selection) were created in this way, and relationships among sires were included in the analyses of these data. In a few cases, the simulated data comprised unrelated base population sire families only, and the data were analysed accordingly.
The data at the phenotypic scale were generated as follows: where the value of the known threshold t, was set depending on the desired level of incidence.
The simulations of the original data were carried out such that the resulting structures led to cases with different amounts of information associated with the parameters of the model. The simulation results were all based on 50 replicates.

RESULTS
Estimates of heritability using MML and GSR (mean of the marginal posterior distribution of heritability) for three data structures are shown in table I. The analyses based on GSR assumed the proper uniform priors for the 'fixed' effect and the variance component, as described in Model for the analysis of binary responses. The structure Al is one in which there is adequate information in the data to make inferences about all the parameters of the model. There is no ECP in this structure and the incidence is set at 60%. Both MML and GSR yield estimates that on average are in good agreement with the true value. In structure A2 there is adequate information related to fixed effects, but there is relatively less information associated with the random effects. In A2 there are no ECPs. In structure A3 there are ECPs, but there is a large amount of data associated with the non-ECP fixed effects and there is little information about the random effects. As expected, the MML yields estimates of heritability that are biased downwards. GSR does not show signs of bias. Structures A4 and A5 represent cases where there is little information about fixed effects. The percentage of ECPs is higher under structure A5 (65% in A5 versus 25% in A4). Both methods produce positively biased results; the bias using GSR is larger than that using MML. Table II shows estimates of the mean of the marginal posterior distribution of the heritability using GSR under a variety of data structures, and performing the Bayesian analysis on the basis of proper and improper posterior distributions. The point of the results in the table is to draw attention to the type of data structures that can cause GSR to yield biased inferences about the variance component and to illustrate that the bias is not necessarily related to the impropriety of the posterior distribution. In cases Bl and B3, the data comprise unrelated sire families and an improper uniform prior was assigned to the fixed effects and to the variance component in such a way that the conditions required to guarantee that the posterior distribution is proper (Natarajan and McCullogh, 1995) are not met.
Although the ECP is present at a level of 27% in Bl, with a large number of families and a large number of observations per family and per fixed effect, the estimate of the mean of the marginal posterior distribution of heritability over the 50 Monte Carlo replicates shows good agreement with the true value. This result is disturbing since the sampler retrieved a meaningful estimate from a posterior distribution which is not guaranteed to be proper! This point is discussed at length in Casella (1996).
In cases B2 and B3 there is a limited amount of information per random effect and per fixed effect. The ECP is not present. The difference between these two cases is that under B2, GSR was implemented with proper priors for all the parameters, which leads to a proper joint posterior distribution. However, the estimates using GSR are biased in both cases and the size of the bias is similar. Thus the presence of bias does not seem to be necessarily related to the fact that the posterior distribution is improper.
In B4, the analysis was carried out assigning proper priors to all the parameters and there are no ECPs (the marginal posterior distribution of heritability is guaranteed to be proper). In common with cases B2 and B3, there are relatively few observations per fixed and random effect. Relative to cases B2 and B3, case B4 comprises a large number of sire families (4 000). Indeed, the mean square error of the heritability is almost three times smaller than in cases B2 and B3. Here, the marginal posterior distribution of heritability is highly peaked and symmetric, but is shifted to the right relative to the true value of 0.5. One of the 50 replicates was arbitrarily chosen and this distribution is shown in figure 1. Despite the large number of sire families, the absence of ECPs and the fact that the marginal posterior distribution is proper, the bias using GSR persists. Table III shows comparisons of GSR and MJP-IBC. Method GSR was implemented assigning normally distributed priors for the fixed effects and the proper uniform distribution in the interval !0, 1/3! for the prior sire variance.
Cases Cl and C2 have the same data structure as B2 in table II, where there are relatively few observations per 'fixed' effect but no ECPs. Using GSR and assigning a normal prior to the 'fixed' effects is effective in removing the bias. Data structures C3 and C4 are identical to A4 and A5. These are characterised by few observations per 'fixed' effect, but with ECP. Again, assigning a normal prior to the 'fixed' effects produces estimates of heritability in better agreement with the true value. In all these four cases, MJP-IBC yields estimates of heritability in good agreement with the true value. The mean square errors of estimates based on GSR and on MJP-IBC are of similar magnitude.
In data structure C5, the frequency of ECP is 80% and the incidence is 90%. The strategy of using GSR treating fixed effects as random yields estimates that are severely biased upwards. MJP-IBC yields estimates that do not deviate significantly from the true value, with a mean square error two and a half times smaller.

DISCUSSION
The purpose of this work has been to draw attention to the poor frequentist properties of the estimator of the mean of the marginal posterior distribution of variance components in mixed threshold models when there is little information associated with the fixed effects. The simulation results indicate that this bias is not necessarily related to the impropriety of posterior distributions. Indeed, we have shown that biased inferences result in cases when the posterior is guaranteed to be proper, and agreement with the true value was obtained when the propriety of the posterior distributions cannot be guaranteed. Further, this is seen as a problem, because the bias persists in situations when this marginal posterior distribution is highly peaked, symmetric and guaranteed to be proper. In other words, this is not the bias to be expected from a non-linear estimator in a small sample setting.
The type of data structures that can lead to this problem is ubiquitous in animal breeding. This is characterised by models with a large number of fixed effects, with relatively few observations per fixed effect. As we have shown, the biased inferences associated with the variance component are not necessarily related to the extreme category problem, although this can be conducive to then.
The inability to obtain a closed form solution to the marginal posterior distribution of the variance component has precluded us from fully understanding the cause of the bias. The lack of asymptotic independence between the estimator of the fixed effect and of the variance component, either in the context of the generalized estimating equation studied by Zeger et al (1988) or in the context of the likelihood of error contrasts reported by Drum and McCullagh (1993), does not necessarily cast light on the present problem. In a Bayesian setting, one is marginalising with respect to the fixed effects, and provided that one does not stumble into problems of impropriety of distributions, it is not clear to us how assertions that are valid from a likelihood viewpoint can help explain the present results.
A strategy that can be interesting to explore is to run the Gibbs sampler after having integrated the fixed effects out analytically from the joint posterior distribution. The relevance of this is based on the observation that running the sampler conditioning on the true values of the fixed effects, yields correct inferences about heritability (Moreno, unpublished).
Using GSR but treating fixed effects as random, removes much of the bias in inference about variance components, as shown by Hoeschele and Tier (1995).
However, we have illustrated that this alternative does not always work. In such cases, it is useful for the researcher to have a choice of other methods to solve the problem. The iterative bootstrap bias correction is an alternative, computationally simple method, that shows good frequentist behaviour. In this work, we chose to use it in conjunction with the mode of the joint posterior distribution, due to its computational simplicity. We stress, however, that the iterative bias correction can in principle be applied to any method, but if it is used, for example, together with the Gibbs sampler, computing time can become a limiting factor. An important point to stress is that the price to pay for unbiasedness is that estimates can fall outside the parameter space. This will never be the case using GSR treating fixed effects as random, or using approximations to the marginal posterior distribution of the variance component. It would certainly be interesting to extend the study of Tempelman (1994) and to explore the behaviour of alternative computational forms of the Laplacian approximation further.
The present simulation study was based on the sire model where each sire effect is represented by several observations. A more relevant model in an animal breeding context is the individual animal model, where each additive genetic value is associated at the most with one observation. In the case of the animal model, the bias problems discussed in this paper are accentuated, as illustrated in the simulation study of Hoeschele and Tier (1995). The study of Tempelman (1994) also showed poor behaviour of the Laplacian approximation in the case of the animal model. It is precisely in such a situation that methods which control this bias such as the MJP-IBC implemented here may be useful. mation of variance components in threshold mixed models. In: Proc of the 5th World Cong Genet Appl to Livestock Prod, University of Guelph, Guelph,22,[15][16][17][18] Wright S (1934) An analysis of variability in number of digits in an inbred strain of guinea pigs. Genetics 19,  Zeger SL, Liang K, Albert PS (1988) Models for longitudinal data: a generalized estimating equation approach. Biometrics 44, [1049][1050][1051][1052][1053][1054][1055][1056][1057][1058][1059][1060] Zeger SL, Karim MR (1991) Generalized linear models with random effects: a Gibbs sampling approach. J Am Stat Assoc 86, 79-86 APPENDIX This appendix shows the log-conditional posterior distributional forms for the fixed and random effects and the variance component of the binary threshold model. The adaptive rejection sampling algorithm was implemented using these distributional forms.
The log-fully conditional posterior distributional form for the xth fixed effect is: where ci = ! _ 1 if i (j) = x and 0 represents the rest of parameters. 0 0 if i(i) 7! x Log-conditional distributional form for the xth random effect is: Log-conditional distributional form for the variance component is: which is proportional to a scaled inverted chi-square distribution.