Inference about multiplicative heteroskedastic components of variance in a mixed linear Gaussian model with an application to beef cattle breeding

Summary - A statistical method for identifying meaningful sources of heterogeneity of residual and genetic variances in mixed linear Gaussian models is presented. The method is based on a structural linear model for log variances. Inference about dispersion parameters is based on the marginal likelihood after integrating out location parameters. A likelihood ratio test using the marginal likelihood is also proposed to test for hypotheses about sources of variation involved. A Bayesian extension of the estimation procedure of the dispersion parameters is presented which consists of determining the mode of their marginal posterior distribution using log inverted chi-square or Gaussian distributions as priors. Procedures presented in the paper are illustrated with the analysis of muscle development scores at weaning of 8575 progeny of 142 sires in the Maine-Anjou breed. In this analysis, heteroskedasticity is found, both for the sire and residual components of variance.

(Received 28 April 1992 ; accepted 23 September 1992) Summary -A statistical method for identifying meaningful sources of heterogeneity of residual and genetic variances in mixed linear Gaussian models is presented. The method is based on a structural linear model for log variances. Inference about dispersion parameters is based on the marginal likelihood after integrating out location parameters. A likelihood ratio test using the marginal likelihood is also proposed to test for hypotheses about sources of variation involved. A Bayesian extension of the estimation procedure of the dispersion parameters is presented which consists of determining the mode of their marginal posterior distribution using log inverted chi-square or Gaussian distributions as priors. Procedures presented in the paper are illustrated with the analysis of muscle development scores at weaning of 8575 progeny of 142 sires in the Maine-Anjou breed. In this analysis, heteroskedasticity is found, both for the sire and residual components of variance. heteroskedasticity / mixed linear model / Bayesian technique R.ésumé -Inférence sur une hétérogénéité multiplicative des composantes de la variance dans un modèle linéaire mixte gaussien: application à la sélection des bovins à viande. Une méthode statistique est présentée, capable d'identifier les sources significatives d'hétérogénéité de variances résiduelles et génétiques dans un modèle linéaire mixte gaussien. La méthode est fondée sur un modèle structurel de décomposition du logarithme des variances. L'inférence concernant les paramètres de dispersion est basée sur la vraisemblance marginale obtenue après intégration des paramètres de position. Un

INTRODUCTION
One of the main concerns of quantitative geneticists lies in evaluation of individuals for selection. The statistical framework to achieve that is nowadays the mixed linear model (Searle, 1971), usually under the assumptions of normality and homogeneity of variances. The estimation of the location parameters is performed with BLUE-BLUP (Best Linear Unbiased Estimation-Prediction), leading to the well-known Mixed Model Equations (MME) of Henderson (1973), and REML (acronym for REstricted -or REsidual-Maximum Likelihood) turns out to be the method of choice for estimating variance components (Patterson and Thompson, 1971): However, heterogeneous variances are often encountered in practice, eg for milk yield in cattle (Hill et al, 1983;Meinert et al, 1988;Dong and Mao, 1990;Visscher et al, 1991;Weigel, 1992) for meat traits in swine (Tholen, 1990) and for growth performance in beef cattle . This heterogeneity of variances, also called heteroskedasticity (McCullogh, 1985), can be due to many factors, eg management level, genotype x environment interactions, segregating major genes, preferential treatments (Visscher et al, 1991).
Ignoring heterogeneity of variance may reduce the reliability of ranking and selection procedures although, in cattle for instance, dam evaluation is likely to be more affected than sire evaluation (Hill, 1984;Vinson, 1987;Winkelman and Schaeffer, 1988).
To overcome this problem, 3 main alternatives are possible. First, a transformation of data can be performed in order to match the usual assumption of homogeneity of variance. A log transformation was proposed by several authors in quantitative genetics (see eg Everett and Keown, 1984;De Veer and Van Vleck, 1987;Short et al, 1990, for milk production traits in cattle). However, while genetic variances tend to stabilize, residual variances of log-transformed records are larger in herds with the lowest production level (De Veer and Van Vleck, 1987;Boldman and Freeman, 1990;Visscher et al, 1991). ' The second alternative is to develop robust methods which are insensitive to moderate heteroskedasticity (Brown, 1982).
The last choice is to take heteroskedasticity into account. Factors (eg region, herd, year, parity, sex) to adjust for heterogeneous variances can be identified. But such a stratification generates a very large number of cells (800 000 levels of herd x year in the French Holstein file) with obvious problems of estimability. Hence, it is logical to handle unequal variances in the same way as unequal means, ie via a modelling (or structural) approach so as to reduce the parameter space, by appropriate identification and testing of meaningful sources of variation of such variances.
The model for the variance components is described in the Model section. Model fitting and estimation of parameters based on marginal likelihood procedures are presented in the Estimation of Parameters, followed by a test statistic in Hypothesis Testing. A Bayesian alternative to maximum marginal likelihood estimation is presented in A Bayesian Approach to a Mixed Model Structure In the Numerical application section, data on French beef cattle are analyzed to illustrate the procedures given in the paper. Finally, some comments on the methodology are made in the Discussion and Conclusion. Following Foulley et al (1990  , the population is assumed to be stratified into I subpopulations, or strata (indexed by i = 1, 2, ... , I) with an (n i x 1) data vector y i , sampled from a normal distribution having mean i ii and variance R. i = a2 ei I&dquo;

MODEL
i . Given ii i and R i Following Henderson (1973), the vector II i is decomposed according to a linear mixed model structure: where X i and Z; are (n i x p) and (n i x q i ) incidence matrices, corresponding to fixed J3 (p x 1 ) and random u i (q i x 1 ) effects respectively. Fixed effects can be factors or covariates, but it is assumed in the following that, without loss of generality, they represent factors.
In the animal breeding context, u i is the vector of genetic merits pertaining to breeding individuals used (sires spread by artificial insemination) or present (males and females) in stratum i. These individuals are related via the so-called numerator relationship matrix A i , which is assumed known and positive definite (of rank q i ).
Elements of u i are not usually the same from one stratum to another. A borderline case is the &dquo;animal&dquo; model ((auaas and Pollak, 1980) where animals with records are completely different from one herd to another. Nevertheless, such individuals are genetically related across herds. Therefore, model [3] has to be refined to take into account covariances among elements of different u!s. As proposed by Gianola et al (1992), this can be accomplished by relating Ui to a general q x 1 vector u * of standardized genetic merits, via the q i x q S i matrix: with A being the overall relationship matrix of rank q, relating the q breeding I animals involved in the whole population, with q x L q j .
i=l Thus, S i is an incidence matrix with 0 and 1 elements relating the q i levels of u * present in the ith subpopulation to the whole vector (q x 1) of u elements. For instance, if stratification is made by herd level, the matrices S i and S i' (i ! i') do not share any non-zero elements in their columns, since animals usually have records only in one herd. On the contrary, in a sire model, a given sire k may have progeny in 2 different herds (i, i') thus resulting in ones in both kth columns of S i and Si.
Notice that in this model, any genotype x stratum interaction is due entirely to scaling  The approach taken here comes from the theory of generalized linear models involving the use of a link function so as to express the transformed parameters with a linear predictor (McCullagh and Nelder, 1989). For variances, a common and convenient choice is the log link function (Aitkin, 1987;Box and Meyer, 1986;Leonard, 1975;Nair and Pregibon, 1988): where wey and w' . are incidence row vectors of size k e and k u , respectively, corresponding to dispersion parameters f g and !u. These incidence vectors can be a subset of the factors for the mean in (2!, but exogeneous information is also allowed. Equations [6] and [7] define the variance component models.
These models can be rewritten in a more compact form as follows.
The vector 0 includes p + q location parameters. The matrix T can be viewed as an &dquo;incidence&dquo; matrix, but which depends here on the dispersion parameters T u through the variances Q ua.
Both variance models can also be compactly written as: The k e + k u dispersion parameters !e and y! can be concatenated into a vector ( T = (T!, T!)' with corresponding incidence matrix W = W e EÐ W u' The dispersion model then reduces to: where a 2 = (CF e 2&dquo; cF! 2' )' and 1n a 2 is a symbolic notation for (In a; 1 ' ... Inaejl 2 In a!1 ' , .. , In a![)'.

ESTIMATION OF PARAMETERS
In sampling theory, a way to eliminate nuisance parameters is to use the marginal likelihood (Kalbfleisch, 1986). &dquo;Roughly speaking, the suggestion is to break the data in two parts, one, part whose distribution depends only on the parameter of interest, and another part whose distribution may well depend on the parameter of interest but which will, in addition, depend on the nuisance parameter. !...! This second part will, in general, contain information about the parameter of interest, but in such a way that this information is inextricably mixed up with the nuisance parameter&dquo; (Barnard, 1970). Patterson and Thompson (1971) used this approach for estimating variance components in mixed linear Gaussian models.
Their derivations were based on error contrasts. The corresponding estimator (the so-called REML) takes into account the loss in degrees of freedom due to the estimation of location parameters.
Alternatively, Harville (1974) proved that REML can be obtained using the noninformative Bayesian paradigm. According to the definition of marginalization in Bayesian inference (Box and Tiao, 1973;Robert, 1992), nuisance parameters are eliminated by integrating them out of the joint posterior density.
Keeping in mind that the sampling and the non-informative Bayesian approaches give rise to the same estimation equations, we have chosen the Bayesian techniques for reasons of coherence and simplicity. The parameters of interest are here the dispersion parameters r, and the location parameters 6 appear to be nuisance parameters. Inference is hence based on the log marginal likelihood L( T ; y) of r: An estimator y of T is given by the mode of L( T ; y): where r is a compact part of R k e +k u.
This maximization can be performed using a result by Foulley et al (1990,1992) which avoids the integration in [13]. Details can be found in the A PP endix. This procedure results in an iterative algorithm. Numerically, let [t] denote the iteration t; the current estimate 9 [Hl] of r is computed from the following system: where i lt] is the current estimate at iteration t, W the incidence matrix defined in !12!, Q M is the weight matrix depending on 0 and on ê [t] , which are the solution and the inverse coefficient matrix respectively of the current system in 0 (this system is described next), z! is the score vector depending on 6 and C!.
Elements of Q l ' l and i lt ) are given in the Appendix. Regarding computations involved in !15!, 2 types of algorithms can be considered as in San Cristobal (1992). A second order algorithm (Newton-Raphson type) converges rapidly and gives estimates of standard errors of y, but computing time can be excessive with the large data sets typical of animal breeding problems. As shown in Foulley et al (1990), a first order algorithm can be easily obtained by approximating the (a matrix in [15] by its expectation component (Qa!,E in the appendix notations). This EM (Expectation-Maximization; Dempster et al, 1977) algorithm converges more slowly, but needs fewer calculations at each iteration and, on the whole, less total CPU time for large data sets.

HYPOTHESIS TESTING
An adequate modelling of heteroskedasticity in variance components requires a procedure for hypothesis testing. Let  to dispersion parameters (ln a 2 = W7 ) ie proceed as if T had a mixed model structure (Garrick and Van Vleck, 1987). To overcome the difficulty of a realistic interpretation of fixed and random effects for conceptual populations of variances from a frequentist (sampling) perspective, one can alternatively use Bayesian procedures. It is then necessary to place suitable prior distributions on dispersion parameters and follow an informative Bayesian approach.
In linear Gaussian methodology, theoretical considerations regarding conjugate priors or fiducial arguments lead to the use of the inverted gamma distribution as a prior for a variance a 2 (Cox and Hinkley, 1974;Robert, 1992). Such a density depends on hyperparameters 77 and s 2 . The former conveys the so-called degrees of belief, and the latter is a location parameter. The ideas briefly exposed in the following are similar to those described in Foulley et al (1992).
Hence, a prior density for y = ln Q 2 can be obtained as a log inverted gamma density. As a matter of fact, it is more interesting to consider the prior distribution of v = &dquo;y &mdash; T°, with q° = In s 2 , ie where r(.) refers to the gamma function. Let us consider a K-dimensional &dquo;random&dquo; factor v such that Vk 1 77k (k = 1, ... K) is distributed as a log inverted gamma InG-l ( 1]k ) ' Since the levels of each random factor are usually exchangeable, it is assumed that 1]k = 1] for every k in {1, ... K}: For v k in [20] small enough, the kernel of the product of independent distributions having densities as in [19] can be approximated (using a Taylor expansion of [19] about v equal to 0) by a Gaussian kernel, leading to the following prior for v: As explained by Foulley et al (1992), this parametrization allows expression of the T vector of dispersion parameters under a mixed model type form. Briefly, from [19] one has 1 = 1° + v or 1 = P 'oS + v if one writes the location parameter -to = In S 2 as a linear function of some vector 8 of explanatory variables (p' being a row incidence vector of coefficients). Extending this writing to several classifications in v leads to the following general expression: where P and Q are incidence matrices corresponding to fixed effects E and random effects v, respectively, with [20] or [21] as prior distribution for v.
Regarding dispersion parameters T , it is then possible to proceed as Henderson (1973) did for location parameters 11 , ie describe them with a mixed model structure. Again, as illustrated by formula [22], the statistical treatment of this model can be conveniently implemented via the Bayesian paradigm.  (Berger, 1985;Robert, 1992): Then, the dispersion parameters are obtained by the mode of the posterior density of T given the hyperparameters equal to their estimates: or similarly for t.  (1992). The extension to the u-part is straightforward.

NUMERICAL APPLICATION
Sires of French beef breeds are routinely evaluated for muscular development (MD) based on phenotypic performance of their male and female progeny. Qualified personnel subjectively classify the calves at about 8 months of age, with MD scores ranging from 0 to 100. Variance components and sire genetic values are then estimated by applying classical procedures, ie REML and BLUP (Henderson, 1973;Thompson, 1979), to a mixed model including the random sire effect and a set of fixed effects described in table I. The second factor listed in table I, condition score (&dquo;Condsc&dquo;), accounts for the previous environmental conditions ( eg nutrition via fatness) in which calves have been raised.
Some factors among those described in table I may induce heterogeneous variances. In particular, different classifiers are expected to generate not only different MD means, but different MD variances as well. Thus, the usual sire model with assumption of homogeneous variances may be inadequate. This hypothesis was tested on the Maine-Anjou breed. After elimination of twins and further editing described in table I, the Maine-Anjou file included performance records on 8 575 progeny out of 142 sires (&dquo;Sire&dquo;) recorded in 5 regions (&dquo;Region&dquo;) and 7 years (&dquo;Year&dquo;). Other factors taken into account were: sex of calves (&dquo;Sex&dquo;), age at scoring (&dquo;Age&dquo;), claving parity (&dquo;Parity&dquo;), month of birth (&dquo;Month&dquo;) and classifier ( &dquo;Classi&dquo; ). In most strata defined as combinations of levels of the previous factors, only one observation was present.

Preliminary analysis
A histogram of the MD variable can be found in figure 1. The distribution of MD seems close to normality, with a fair PP-plot (although the use of this procedure is somewhat controversial), and skewness and kurtosis coefficients were estimated as -0.09 and 0.37 respectively. Some commonly used tests for normality rejected the null hypothesis, while others did not reject it, namely Geary's u, Pearson's tests for skewness and kurtosis (Morice, 1972)  The usual sire model with all factors from table I in the mean model, and variance components estimated by EM-REML, was fitted, leading to estimates 6d = 70.1l,a,2, = 6.91, and h 2 = 46fl /(6d + 3!) = 0.36. Note that this model is equivalent, in our notation, to the homogeneous model in f g and Yu .

Search for a model for the variances
The following additive mean model M B was considered as true throughout the whole analysis This model was chosen in agreement with technicians of the Maine-Anjou breed and is used routinely for genetic evaluation of Maine-Anjou sires.
A forward selection of factors strategy was chosen to find a good variance model My but in 2 stages; a backward selection strategy would have been difficult to implement because of the large number of models to compare and the small amount of information in some strata generated by those models.
(i) since a2 represents > 90% of the total variation, it was decided to model that component first, assuming the r u -part homogenous; ' (ii) the &dquo;best&dquo; T u -model was thereafter chosen while keeping unchanged the &dquo;best&dquo; T e-model. The different nested models were fitted using the maximum marginal likelihood ratio test (MLRT) A described in !17J.
During the first stage (i), the homogeneous sire variance was estimated, for computational ease, with an EM-REML algorithm, and the Te parameter estimates were calculated as in Foulley et al (1992). This strategy leads, of course, to the same results as those obtained with the algorithm described in the Estimation of parameters.
The first step consisted of choosing the best one-factor variance model from results presented in table lib. The next steps, ie the choice of an adequate 2-factor model, and then of a 3-factor model, etc, are summarised in table III. Finally, the following additive model was chosen: The model can also be simplified after comparing estimates of factor levels, and then collapsing these levels if there are not significantly different.
For the (ii) stage, the &dquo;best&dquo; r u -model was the model (see table IV): We were not able to reach convergence of the iterative procedure for the models (Mo, M.y e , Classi) and (Mo, M!(,, Region), although some levels of the Classi factor were collapsed. This phenomenon is related to a strong unbalance of the design: for instance, one classifier noted the calves of only 4 sires, making quite impossible a coherent estimation of Classi-heterogeneous sire variances. The other factors (except Year) had no significant effect on the variation of the sire variances. Because of imbalance, the model gave unsatisfactory results eg heritability estimates greater than one. !e,Month = 0.0024 respectively, or alternatively using % values of the coefficient of variation for ae, (!e CV 2 ) CV e , Class i = 14.5%, CV,, Yar = 9.5% and CV e , Month = 4.9% respectively. In fact, the smaller the cell size (n i ), and the smaller CV, the greater the shrinkage of the sample estimates (6 f ) toward the mean variance (3 ) since the regression coefficient toward this mean in the equation Q2 = õ'2 + b(6 i 2 _ õ'2) is approximately b = n d[ n i + (2/CV!)] with !7 = 2/CV 2 : see also Visscher and Hill (1992).

Estimates of the dispersion parameters
The genetic variation in heifers turns out to be less than one half what it is in bulls even though the phenotypic variance was virtually the same. This may be due to the fact that classifiers do not score exactly the same trait in males (muscling) as in females (size and/or fatness). It may also suggest that the regime of male calves is supplemented with concentrate. Location parameters are compared in figures 2a-d under different dispersion models, through scatter plots of estimates of standardized sire merits (u * ). Indexes based on &dquo;subclass means&dquo; (V i = y i , i = 1, ... I, with homogeneous variances) and those based on the &dquo;sire model&dquo; under the homogeneity of variance assumption are far away from each other (see fig 2a). Figure 2a is just a reference of discrepancy, which illustrates the impact of the BLUP methodology. When heterogeneity is introduced among residual variances, sires' genetic values do not vary too much, as shown in figure 2b. Modelling of the genetic variances has a larger impact on the sire genetic values (see figure 2c) than modelling of residual variances. Finally, the Bayesian treatment of r e -parameters by introducing random effects in the model (M B , M&dquo; YJ does not have any influence on the sire genetic merits (fig 3d).
Evaluation of sires can be biased if true heterogeneity of variance is not taken into account. As shown in table VI, sire number 13 went down from the 16th to the 24th position because his calves were scored mostly by classifier no 1 who uses a large scale of notation (see T -estimates in table V). On the other hand, sire 103 went up from the 25th to the 14th place since the corresponding Classi and Condsc levels have low residual variance (for the other factor levels represented, the variances were at the average). For the same reason, the sire genetic merits were also affected by modelling In ad. The difference in genetic merit for sire 56 (1.40 vs 1.74 under the homoskedastic and the residual heteroskedastic models respectively) is also explained by the fact that the calves of this sire were scored exclusively by classifier no 12 and in 1983 (Year = 1). Due to modelling Q u, this sire went down again (from 1.74 to 1.63 under the full heteroskedastic model) because all its progeny are females with a lower Q u component than in males. Other things being equal, a reduction in the oru variance results in a larger ratio, or equivalently a smaller heritability and consequently in a higher shrinkage of the estimated breeding value toward the mean. In other words, if a decrease in genetic variance is ignored, sires above the mean are overevaluated and sires below the mean are underevaluated.

Hypothesis checking
Normality assumptions made in [1] and [5] were checked at each step of the analysis. After modelling residual variances, the distribution of standardized residuals became closer to normality, in terms of skewness and especially kurtosis. This phenomenon was observed in the whole sample and also in the subsamples defined by the levels of the factor considered in r e . On the other hand, normality of the residuals was stable in the subsamples defined by the factors absent from the r emodel.
Normality of the distribution of the standardized sire values in terms of kurtosis and PP-plot was continuously damaged at each step of the variance modelling: estimated kurtosis was 0.61, 0.72 and 0.90, for the homoskedastic, residual heteroskedastic and fully heteroskedastic models respectively. Moreover, skewness for the 142 sire genetic merits improved slightly during that process: -0.09, -0.003 and -0.03 for the same models respectively.

Computational aspects
Programmes were written in Fortran 77 on an IBM 3090 by implementing an EM algorithm corresponding to [15]. The convergence was fast: 15-20 cycles for heteroskedastic T e-models with < 7J estimated by EM-R.EML ((i) stage), and 15-40 cycles for fully heteroskedastic T -models or heteroskedastic T e -models with random effects. CPU time was between 2-5 min per model fit (estimation of parameters and computation of the log marginal likelihood.

DISCUSSION AND CONCLUSION
This paper is an extension to u-components of variances of the approach developed by Foulley et al (1992) to consider heterogeneity in residual variances using a structural model to describe dispersion parameters, in a similar way as usually done on subclass means.
In that respect, our main concern focuses on ways to render models as parsimonious as possible so as to reduce the number of parameters needed to assess heteroskedasticity of variances. An interesting feature of this procedure is to assess, through a kind of analysis of variance, the effects of factors marginally or jointly. For instance, one can test heterogeneity of sire variances among breeds of dams after adjusting for possible sources of variation such as management level. In the same way, differences among group of sires in within-sire variances (which might be related to a segregating major gene) can be tested while taking into account the influence of other nuisance factors (season, nutrition...). However, the power of the likelihood ratio test for detecting heterogeneity of variance can be a real issue in many practical instances.
From the genetic point of view, the approach is quite general since it can deal with heterogeneity among within and between family components of variances, or among genetic and environmental variances. Factors involved for u and e components of variance may be different or the same, making the method especially flexible. Our modelling allows one to assume (or even test) whether the ratios of variances or heritabilities are constant over levels of some single factor or combination of factors (Visscher and Hill, 1992). If a constant heritability or ratio of variances a = or 2 i/ or 2 among strata is assumed, the model involves the parameters y and a only, and reduces to In o, ei 2 = we!re with oru 2i replaced by a; j 0: in the likelihood function.
The shrinkage estimator for the variances proposed by eg, Gianola et al (1992), follows the same idea of the Bayesian estimator described in the Bayesian approach section. When a Gaussian prior density is employed for the dispersion parameters Y , the hyperparameter 6 acts as a shrinker. But the Bayesian approach for a direct shrinkage of variance components assumes that heterogeneity in such components (residual and u components) is due only to one factor. The approach presented in this paper is more general since it can cope with more complex structures of stratification which may differ from one component to the other. Moreover, its mixed model structure allows great flexibility to adjust variances in relation to the amount of information for factors in the model; eg when data provide little information for some factors (or levels) or considerable for others, our procedure behaves like BLUP (or James-Stein) ie shrink estimates of dispersion parameters toward zero if there is little information; only with sufficient information can the estimate deviate. For instance, our methodology provides a simple and rational procedure to shrink herd variances (whatever they are, genetic, residual or phenotypic) towards different population values (eg regions, as proposed by Wiggans and VanRaden, 1991) due to poor accuracy of within herd or herd-year variances (Brotherstone and Hill, 1986). It then suffices to use a hierarchical (linear) mixed model for herd log-variances and take the population factor ( eg region) as fixed and herd as random within that factor. An illustration of the flexibility and feasibility of our procedure was recently given by Weigel (1992) in analyzing sources of heterogeneous variances for milk and fat yield in US Holsteins.
Coming back to the case of a unique factor of variation for the sire variances, Although this presentation is restricted to a single random factor u * , it can be generalized to a multiple random factor situation. If such factors are uncorrelated, the extension is straightforward. When covariances exist, one may simply assume, as proposed by Quaas et al (1989), that heterogeneity in covariances is due to scaling. This means, for instance, that in a sire (s i )-maternal grand sire (t!) model Y ijk = X!. k 13 + Si +t! +e2!k, one will model as h' 2 a th 2 as previously, and assume that the covariance is Q9t ,, = pa sha t h for stratum h. If the model is parameterized in terms of direct a o and maternal a m effects as follows through the transformation s i = ao i/ 2 and t j = ao!/4 + a m; /2, one can set the genetic correlation p a to a constant, ie a ao = P a a aoh O a , n -Notice that this condition is not equivalent to the previous one, except if a aoh /!d&dquo;,h does not depend on h.
Although the methodology is appealing, attention must be drawn to the feasibility of the method. The first problem is the inversion of the coefficient matrix in [16] required for the computation of the variance system (15]. In animal breeding applications, this matrix is usually very large. This limiting factor is already becoming less important due to constant progress in computing software and hardware. The technique of absorption is usually used to reduce the size of matrices to invert. Another approach is to approximate the inverse. One can, for instance, use a Taylor series expansion of order N for a square invertible matrix A where the square matrix A o is a matrix close to A and is, of course, easy to invert, and where 11 -11 denotes some norm on the space of invertible matrices. Methods viewed in Boichard et al (1992) can also help to approximate A-1 in particular cases such as sparse matrices, &dquo;animal model&dquo;, etc.
Statistical power for likelihood ratio tests was investigated for detection of heterogeneous variances in the usual designs of quantitative genetics and animal breeding. Results given by Visscher (1992) and Shaw (1991) indicate generally low power values for detecting heterogeneity in genetic variance. According to Shaw, a nested design of 900 individuals out of 100 sire families provides a power of 0.5 for genetic variances differing by a factor of 2.5. This clearly indicates the minimum requirements in sample size and family numbers which should be met before carrying out such an analysis and the limits therein. Therefore it seems unrealistic to model genetic variances in practice according to more than 1 or 2 factors, and it might be wise to consider some of them as random if little information is provided by the data in each level of such factors.
Although  (Box and Tiao, 1973). This problem does not occur with an &dquo;animal model&dquo;, but can arise when a &dquo;sire model&dquo; is used, and is not specifically related to heteroskedasticity.
From a statistical point of view, the procedure uses the concept of variance function (Davidian and Carroll, 1987) as an extension to dispersion parameters of the link function. Our presentation focuses on the log link function which is the most common choice in this field (see for instance San Cristobal, 1992, for a review of variance models) &dquo;for physical and numerical reasons&dquo; (Nair and Pregibon, 1988). Following Davidian and Carroll (1987) or Duby et al (1975), the question can be asked whether or not variances vary according to means or location parameters. In the Maine-Anjou data, however, it does not seem to be the case, thus validating our choice in !10!.
It would be interesting to extend our method to a fully generalized linear mixed model on means and on variances with or without common parameters between the mean model and the variance model. Numerical integration or Gibbs sampling procedures would then be required although approximate methods of inference can also be used for such models (Breslow and Clayton, 1992;Firth, 1992). Statistical problems arising with common parameters are already highlighted by van Houwelingen (1988).
With a fully fixed effect variance model, techniques of estimation and hypothesis testing for dispersion parameters presented here are those of the classical theory of likelihood inference (likelihood and likelihood ratio test), except that the marginal likelihood function L( Y ;y) was preferred to the usual likelihood L(13, T ; y), in the light of ideas behind REML estimators of variance components. This test reduces to Bartlett's test (Bartlett, 1937) for a one classification model in variances and under a saturated fixed model on the means (ie ji i = y i , i = 1, ... I). Unfortunately, Bartlett's test is known to be sensitive to departure from normality (Box, 1953).
Simulations are needed to study the robustness of this test and other competing tests. From a Bayesian perspective, the Bayes factor is usually applied for hypothesis testing (see Robert, 1992, for a discussion). The posterior Bayes factor (Aitkin, 1991) could also be used to compare dispersion models, but numerical integration would then be required (see the expression of the likelihood in !18!).
In this paper, focus was on an appropriate way to model heterogeneous variances, but the initial motivation was a best fitting of location parameters (animal evaluation for animal breeders). This difficult problem of feedback, also related to the Behrens-Fisher problem, has to be solved in our particular approach.
Moreover, a great research perspective is open on the important and complicated question of the joint modelling of means and variances (Aitkin, 1987;Nelder, 1991;Helder and Lee, 1991).