Analysis of a genetically structured variance heterogeneity model using the Box–Cox transformation

Summary

Over recent years, statistical support for the presence of genetic factors operating at the level of the environmental variance has come from fitting a genetically structured heterogeneous variance model to field or experimental data in various species. Misleading results may arise due to skewness of the marginal distribution of the data. To investigate how the scale of measurement affects inferences, the genetically structured heterogeneous variance model is extended to accommodate the family of Box–Cox transformations. Litter size data in rabbits and pigs that had previously been analysed in the untransformed scale were reanalysed in a scale equal to the mode of the marginal posterior distribution of the Box–Cox parameter. In the rabbit data, the statistical evidence for a genetic component at the level of the environmental variance is considerably weaker than that resulting from an analysis in the original metric. In the pig data, the statistical evidence is stronger, but the coefficient of correlation between additive genetic effects affecting mean and variance changes sign, compared to the results in the untransformed scale. The study confirms that inferences on variances can be strongly affected by the presence of asymmetry in the distribution of data. We recommend that to avoid one important source of spurious inferences, future work seeking support for a genetic component acting on environmental variation using a parametric approach based on normality assumptions confirms that these are met.

1. Introduction

The classical model of quantitative genetics assumes that genotypes differ in mean but the environmental variance (conditional variance of phenotype, given genotype) is homogeneous across genotypes (this model is referred to as HOM hereinafter). The assumption of homogeneity may not be valid and models with non-genetically structured departures from variance homogeneity were introduced in the 1990s (Foulley et al., 1992; Gianola et al., 1992; San Cristobal et al., 1993). In particular, Foulley & Quaas (1995) proposed models of heterogeneity for both residual and other components of variance.

A significant extension of the classical genetic model, with implications in breeding and evolutionary studies, posits that both mean and variability differ between genotypes (San Cristobal-Gaudy et al., 1998). (The extended model is labelled HET hereinafter). In animal and plant breeding, the focus is to change the mean of a trait by selection; the HET model opens the possibility to reduce variation by selection leading to more homogeneous products (e.g. Hill & Zhang, 2004; Mulder et al., 2007, 2008). In evolutionary biology, most of the models developed to study the maintenance of phenotypic variation assume that environmental variance is constant and explain the levels of variation by invoking a balance between the gain of genetic variance by mutation and its loss by different forms of selection and drift. Zhang & Hill (2005) discuss models where environmental variance is partly under genetic control and study conditions for its maintenance under stabilizing selection.

Early evidence for a genetic component affecting environmental variation stems from comparisons of levels of variation between inbred lines and the F₁ cross between them, with inbreds showing in general larger variance (Falconer & Mackay, 1996). More recently, experimental evidence from isogenic chromosome substitution lines in Drosophila was presented by Mackay & Lyman (2005). Using genetic marker information, support for genetic variance heterogeneity has been reported also by Ordas et al., (2008) in maize and by Ansel et al. (2008) in yeast. Ansel et al. (2008) work with isogenic yeast cells and provide convincing evidence for heterogeneity of gene expression between genotypes, and identify three quantitative trait loci (QTLs) involved in its control.

With the exception of the limited experimental evidence specifically designed to test for genetic control of residual variability, support for the presence of genes acting on environmental variation can be sought fitting the HET model to data and studying the quality of its fit using modern computation tools. Evidence of this kind has been reported using litter size in pigs (Sorensen & Waagepetersen, 2003), adult weight in snails (Ros et al., 2004), uterine capacity in rabbits (Ibáñez et al., 2008), body weight in poultry (Rowe et al., 2006; Wolc et al., 2009), slaughter weight in pigs (Ibáñez et al., 2007) and litter size and weight at birth in mice (Gutierrez et al., 2006).

In all the cases where the quality of fit of the HET and HOM models was compared, various statistical criteria favoured the HET model. An important issue is to clarify whether these results are an artefact of the scale of measurement. Under the HET model, the skewness of the marginal distribution of the data is directly proportional to the correlation between additive genetic values affecting mean and variance (Ros et al., 2004). However skewness can arise in other ways and if it does it can cause the HET model to fit better than the HOM model, leading to the spurious conclusion that genetic heterogeneity is present, when it is absent. Alternatively, genetic heterogeneity of variance may be present, in which case skewness from other causes may affect inferences about the correlation between genetic effects on mean and variance.

These problems can be partly addressed finding the appropriate scale on which to analyse the data using the Box–Cox transformation parameter (Box & Cox, 1964). A large body of literature dealing with transformations has developed following the classical paper of Box & Cox (1964). Early work using transformations in animal breeding are Ibe & Hill (1988) and Gianola et al. (1990). Compelling reasons for using the transformed model have typically been the need to achieve additivity of parameters at the level of the mean, constancy and homogeneity of the error variance and normality of the distribution of observations, although a transformation may not achieve all these goals simultaneously. Here, a transformation is sought that induces normality and linearity at the level of the conditional distribution of data. This is an important modelling assumption that must be fulfilled to avoid distorted inferences.

The main objective of this work is to investigate whether an analysis performed under the ‘appropriate’ scale still supports the HET model, or whether simpler models are to be preferred. The Box–Cox transformation parameter is regarded as an unknown to be inferred from the data, jointly, with the remaining parameters of the model. The model is fitted to data on litter size in pigs and rabbits. In order to obtain insight into the quality of inferences about the transformation parameter when it is fitted simultaneously with the remaining parameters, including the coefficient of correlation, a simulation study involving sets of data that differ in the number of repeated records per individual is also conducted.

Implementation of the model is via McMC (Markov chain Monte Carlo) and two steps are involved. The first step aims at finding the modal value of the marginal posterior distribution of the transformation parameter. In the second step, the model is fitted conditional on this modal value. Several McMC strategies were investigated and the reported results are based on a fixed scan hybrid algorithm that resulted in the best mixing behaviour.

The article is organized as follows: section 2 introduces the Box–Cox model with genetically structured variance heterogeneity and discusses the prior specifications and McMC implementation. Section 3 describes briefly the simulated, the rabbit and the pig litter size data. Section 4 contains the outcome of applying the methodology to the simulated data and to the rabbit and pig litter size data. A discussion is provided in section 5.

2. Methods

(i) Box–Cox model with genetically structured variance heterogeneity

The proposed model, labelled BCHET hereinafter, combines the HET model with a Box–Cox transformation parameter. Let y=(y _i)_i=1ⁿ be the vector of records, which are assumed to be strictly positive. The model states that conditionally on the vectors of location and dispersion parameters, for some unknown transformation parameter λ, the distribution of the transformed vector of phenotypic data y^(λ)=(y _i^(λ))_i=1ⁿ satisfies the normal theory assumptions, that is,

(1)

Here R, the conditional variance of the sampling model, is the diagonal matrix with on its diagonal, and is interpreted as the environmental variance. The transformation is defined as follows (Box & Cox, 1964):

(2)

Note that here the transformation parameter λ is introduced to guarantee normality of the conditional distribution and linearity of the conditional expectation, but not to achieve constant variance. Indeed, R in (1) is heterogeneous due to the presence of , a* and p*; the last two terms being specific to each individual.

Matrices X, , Z and W are known incidence matrices, the scalar μ is the mean, the vector b contains the regression coefficients, the scalar exp (μ*) is the mean environmental variance, while the vector exp (b*) contains the environment-specific scaling effects. The vectors of genetic effects a=(a_i )_i=1^{n ₁} and are assumed to jointly follow a multivariate normal distribution

where A is the known additive relationship matrix and G is a 2×2 covariance matrix

(3)

In equation (3), ρ is the coefficient of genetic correlation between a and a*, and σ_a² and are the additive genetic variances at the level of the mean and environmental variance, respectively. Vectors and represent permanent environmental effects (on mean and on variance) and are assumed to be independently normally distributed with

where I is the identity matrix, and σ_p² and σ_p*² are the permanent environmental variances at the level of the mean and environmental variance, respectively.

The coefficient of skewness of the marginal distribution of the datum y_i ^(λ) (marginalized over a, a*, p, p*) under BCHET is (Ros et al., 2004)

indicating that the sign of the coefficient of skewness is governed by ρ.

(ii) Conditional likelihood function

The conditional likelihood is proportional to the probability density function of the data given the parameters, including a, a*, p, p*. Written in terms of the original untransformed observations, from (1), this leads to

(4)

where

(5)

R and y^(λ) have been defined in connection with (1), and the Jacobian of the transformation is

(iii) Prior specifications

The general size and range of the transformed phenotypic data y^(λ) may depend strongly on λ, and so would the interpretation of the model parameters. Therefore, the prior distributions must be carefully chosen to avoid meaningless inferences. Here, we follow the approach suggested in Box & Cox (1964) to establish the conditional prior distribution, given λ, whereby the relationship between parameters is consistent under different values of λ. Calculations sketched out in Appendix A show that the prior distribution of the parameters that depend on λ is

(6)

(iv) Posterior distribution

Including in (6) the prior distribution of the parameters conditionally independent of λ, , leads to the joint prior distribution for all model parameters. This takes the form

(7)

The joint posterior distribution is proportional to the product of (4) and the prior (7).

The inferential strategy in this study is based on, firstly, estimating the modal value of λ from its marginal posterior distribution. This is done fitting jointly all the parameters of the BCHET model, including the correlation coefficient ρ. The final inferences are based on a second implementation of the BCHET model, where the value of λ is fixed at its posterior mode.

(v) Model comparison and model checking

Models are compared using the pseudo marginal probability of data. This is a method of comparison (Gelfand, 1996) based on the relative global fits of the models. In addition, the skewness of the conditional distribution of the data under both models is investigated using posterior predictive model checking. This provides insight into the ability of the models to address a specific feature of the data that is of primary interest. In this study, one wishes to investigate whether the transformation succeeds in inducing normality at the level of the conditional distribution of the data, which is a basic premise of the model and a necessary condition for correct inferences. Details are in Appendix B.

(vi) McMC algorithm

The BCHET model (1) is implemented using an McMC algorithm, where the components in the model are updated sequentially. The algorithm is based on a combination of Gibbs updates, updates based on random walk proposals and updates based on Langevin–Hastings proposals. In addition, a reparameterization of the additive genetic values (a, a*) to standard normal variates (γ, γ*) was performed to improve the mixing of the algorithm.

In summary, the BCHET model (1) is implemented using the following McMC algorithm.

• • Metropolis–Hasting steps:

  1. 1. (λ, μ, exp (μ*),b) updated jointly using a random walk proposal for λ and fully conditional proposals for (μ, exp (μ*),b).

  2. 2. exp (b*) and (γ, γ*) updated using Langevin–Hastings proposals.

  3. 3. (p, p*) updated using a fully conditional proposal for p and a Langevin–Hastings proposal for p*.

  4. 4. (σ_a², σ_a*²) updated using random walk proposals (on the log scale).

  5. 5. ρ updated using a random walk proposal.

• • Gibbs Steps: σ_p², σ_p*².

Convergence of the algorithm was checked informally using traceplots. We report confidence intervals of Monte Carlo estimates of various posterior means in order to give an idea of the accuracy of the Monte Carlo computations.

3. Data

(i) Simulated data

It is difficult to gain insight about the degree of identifiability of λ in the marginal likelihood when it is fitted jointly with ρ (together with the remaining parameters of the model), because the marginal likelihood cannot be written in closed form. One can conjecture that the number of repeated records per individual plays an important role. Therefore, a small simulation study was conducted to investigate how the quality of inferences about λ from its marginal posterior distribution is affected with increasing number of records per individual. A family structure based on the pedigree of the rabbit litter size data was simulated, with 1, 2, 3, 5 and 10 records per individual. The true value of λ was set equal to 1, and of ρ equal to −0·74, in all cases. These true values are compared with those inferred from the marginal posterior distribution of these parameters.

(ii) Rabbit litter size data

The data originate from a selection experiment for uterine capacity in rabbits spanning 10 generations, carried out at the Universidad Politécnica de Valencia. Uterine capacity will be referred to as litter size hereinafter. Details of the selection experiments and data can be found in Argente et al. (1997) and Ibáñez et al. (2008). The number of individuals in the pedigree was 1281, where 929 of them had litter size records (a maximum of four records per female), and the total number of records was 2996. Reproduction was organized in discrete generations.

(iii) Pig litter size data

The data originate from a large-scale selection experiment for the total number of piglets born per litter (referred to as litter size hereinafter). Details of the selection experiments and data descriptions are provided in Sorensen et al. (2000). The number of individuals in the pedigree was 6437, where 4078 of them had phenotypic records, and the total number of litter size records was 9778.

4. Results

(i) Simulation study

The results reported for different simulated data sets were obtained from McMC output consisting of between 1 000 000 and 3 000 000 iterations. The total number of stored samples was 10 000 for each simulated data set and these were used to infer the marginal distribution of λ.

Results are shown in Table 1. The 95% posterior intervals of λ always include the true value 1 in all simulated data sets, and inferences are sharper with increasing number of records per individual. The posterior distributions are symmetric (not shown); posterior means and modes agree to two decimal figures.

Table 1. MC estimates of posterior means and 95% posterior interval (95% P.I.) for λ and ρ with increasing number of repeated records per individual

Estimates of the posterior means are clearly subject to sampling variation of the data. This frequentist variation may be characterized by replication of the experiment, but this was judged to be computationally too demanding.

The simulation study showed that the posterior correlation of λ and ρ decreased with increasing number of records per individual (not shown). In the data set mimicking the rabbit structure, this posterior correlation is equal to 0·34. The results indicate that the size of the smallest data set, corresponding to the rabbit experiment, is adequate to arrive at satisfactory inferences about features of the marginal posterior distribution of λ.

(ii) Analyses of litter size in rabbits and pigs

(a) The models

The models for the rabbit and pig data are similar to those in Sorensen & Waagepetersen (2003) and Ibáñez et al. (2008). The scalar μ is the mean and the scalar exp (μ*) is the mean environmental variance. In the rabbit data, vector b contains year–season effects with 30 levels and parity order effects with four levels (first, second, third and fourth or higher parities). The vector exp (b*) contains the environment-specific scalings effects associated with year–season and parity. Vectors a and a* contain additive genetic values with 1281 levels, and vectors p and p* contain permanent environmental effects with 929 levels. The parameters of the scaled inverted chi-square distributions associated with σ_a², σ_a*², σ_p², σ_p*², are ν=5, , , and . In the pig data, b is a vector containing the effects of four categorical variables: parity (eight levels), season (four levels), herd (50 levels) and type of insemination (artificial or natural). The vector exp (b*) contains the environment-specific scalings effects with two categorical variables: parity (eight levels) and type of insemination (artificial or natural). Vectors a and a* contain additive genetic effects with 6437 levels, and vectors p and p* contain permanent environmental effects with 4078 levels. The parameters of the scaled inverted chi-square distributions associated with σ_a², σ_a*², σ_p², σ_p*², are ν=5, , , and . The hyperparameters of the prior distributions are taken from Sorensen & Waagepetersen (2003) and Ibáñez et al. (2008).

(b) The mode of λ

Inference about the modal value of λ from its marginal posterior distribution was accomplished running the McMC algorithm with a chain length equal to 1 405 000 in the rabbit data and 750 000 in the pig data. In the rabbit data, the Monte Carlo estimate of the mode of λ was equal to 1·41 and the 95% posterior interval was (1·30, 1·51). In the pig data, the corresponding figures are 1·39 and (1·32, 1·47). In both cases, the posterior distribution was symmetric (the estimated posterior means agreed with the modal values to two decimal figures). In neither case do the 95% posterior intervals of λ contain the value 1. This indicates that in both species the phenotypic data are not conditionally normally distributed under λ=1.

In the rabbit data, the 95% confidence interval of the posterior mean, reflecting Monte Carlo sampling variation, was (1·405, 1·411) and (1·390,1·398) in the pig data, indicating that the degree of accuracy of the Monte Carlo computation was satisfactory.

In order to gain insight into the sensitivity of inferences on prior assumptions, a similar analysis was conducted with scale parameters , , , set equal to 0·124, 0·024, 0·066 and 0·018 in the rabbit data and to 0·234, 0·014, 0·090 and 0·009 in the pig data. Under these prior assumptions, the posterior mode of λ in the rabbit data is 1·40, and the 95% posterior interval for λ was (1·31, 1·51). In the pig data, the corresponding figures are 1·39 and the 95% posterior interval is (1·32, 1·46) reflecting minor differences with the original values in both cases.

(c) Variance components

Results for the rabbit data set are shown in the top rows of Table 2. The figures report Monte Carlo estimates of posterior means, 95% posterior intervals and 95% Monte Carlo confidence intervals for chosen parameters of BCHET under λ=1 (M) and BCHET under the modal value of λ=1·41 (M_λ). Estimates of σ_a² and σ_p² increase by factors of 3 and 7, respectively, when the data are analysed under M_λ, as expected, due to dependence on λ. Despite the lack of prior dependence of σ_a*² and σ_p*² on λ, the posterior mean of σ_a*² decreases by a factor >2 under M_λ, while σ_p*² is slightly reduced. The posterior mean of ρ changes from −0·73 under M to 0·28 under M_λ, and the 95% posterior interval of ρ includes zero under M_λ.

Table 2. MC estimates of posterior means (first row for each model), 95% posterior intervals (second row for each model) and of 95% Monte Carlo confidence intervals (third row for each model) of chosen parameters based on models M and M_λ, for the rabbit litter size data (top) and on models M and M_λ for the pig data set (bottom).

Under M, the smallest posterior mean of the environmental variance (95% posterior intervals in brackets), corresponding to year–season 30 and parity 1, was 3·28 (1·55,5·34). The corresponding largest number, for year–season 15 and parity 3 was 8·23 (4·61, 11·92). Under M_λ, the corresponding figures are 11·16 (5·72, 17·56) and 36·92 (26·39, 52·99).

Under M, there is one heritability for each combination of environment-specific scalings effects at the level of the environmental variance. The average heritability over all combinations of environmental effects is 0·12 with a minimum of 0·09 (0·04, 0·14) (year–season effect 15 and parity 3) and a maximum of 0·18 (0·09, 0·28) (year–season effect 30 and parity 1). Under M_λ, the average heritability is 0·10 with a minimum of 0·07 (0·03, 0·11) (year–season effect 15 and parity 3) and a maximum of 0·19 (0·08, 0·31) (year–season effect 30 and parity 1), reflecting little change due to the transformation. It can be shown that the prior heritability is not affected by the transformation (to the degree of approximation (A1)), a result that agrees with the conclusion from Solomon & Taylor (1999), who worked with a variance homogeneous, one-way classification model.

Figure 1 shows the estimated marginal posterior distributions of σ_a², σ_a*², σ_p² and σ_p*² based on M_λ. The solid line superimposed on each graph represents the densities of the scaled inverted chi-square prior distributions of σ_a², σ_a*², σ_p² and σ_p*². On the top of the figure, the variance parameters σ_a², σ_a*², σ_p² and σ_p*² are assigned scaled inverted chi-squared distributions with degrees of freedom ν=5 and scale parameters , , and , resulting in prior modes equal to 0·35, 0·07, 0·19 and 0·05, respectively. The bottom of the figure shows corresponding distributions obtained setting , , and , leading to prior modes equal to 0·09, 0·02, 0·05 and 0·013, respectively, a four-fold reduction. Prior input does not influence inferences about σ_a² and σ_p² but σ_a*² and σ_p*² are strongly affected. In fact, the posterior distribution of σ_a*² shows little divergence from the prior distribution in both scenarios, suggesting that there is little information in the data about this parameter.

Figure 1. Top: Monte Carlo estimates of σ_a² and σ_{a *}² (left) and σ_p² and σ_p*² (right) under M_λ for rabbit data. The thick lines represent the prior scaled inverted chi-square densities with degrees of freedom ν=5 and scale parameters , , and under M_λ. Bottom: as above, with ν=5 and , , and .

In conclusion, the analysis indicates that the support for a genetically structured variance heterogeneous model in the rabbit data is markedly weaker when the data are analysed under a scale corresponding to the modal value of λ.

Estimated posterior means of dispersion parameters for the pig data set are shown in the bottom rows of Table 2. The figures show a similar pattern as for the rabbit litter size data. The posterior mean of ρ changes sign from negative under M (BCHET with λ=1) to positive under M_λ (BCHET conditional on the mode of λ=1·39). However, in contrast with the rabbit data, the 95% posterior interval of ρ does not include the value ρ=0. Thus, the analysis in the transformed scale results in a strong, positive correlation between additive genetic values affecting mean and variance.

Under M, the smallest posterior mean of the environmental variance, corresponding to natural insemination and parity 1, is 5·87 (5·27, 6·53). The largest posterior mean, corresponding to artificial insemination and parity 2, is 9·33 (8·27, 10·42). Under M_λ, these figures are 33·01 (30·01, 36·02) and 58·06 (47·59, 69·62).

Under M, there is one heritability for each combination of the environment-specific scalings effects . The average heritability is 0·16 with a minimum of 0·14 (0·11, 0·18) (artificial insemination and parity 2) and a maximum of 0·20 (0·16, 0·25) (natural insemination and parity 1). Under M_λ, the average heritability is 0·14 with a minimum of 0·12 (0·08, 0·15) (artificial insemination and parity 6) and a maximum of 0·18 (0·13, 0·23) (natural insemination and parity 1).

Figure 2 illustrates the sensitivity of inferences to prior assumptions for the pig data. In contrast with the rabbit data, a five-fold reduction in prior modes leads to very small changes in the posterior distributions, except for σ_p*².

Figure 2. Top: Monte Carlo estimates of σ_a² and σ_a*² (left) and σ_p² and σ_p*² (right) under M_λ for pig data. The thick lines represent the prior scaled inverted χ² densities with degrees of freedom ν=5 and scale parameters , , and under M_λ. Bottom: as above, with ν=5 and scale parameters , , and .

To give a further idea of the sensitivity of inferences to λ, the rabbit data were also analysed conditional on the lower (λ=1·30) and the upper λ=1·51 cut-off points of the 95% posterior interval for λ. The posterior distribution of σ_a*² was hardly affected, but the posterior distribution of ρ shifted to the right with increasing values of λ. The posterior means and 95% posterior intervals for λ=1·30, 1·41 and 1·51 were −0·10 (−0·55, 0·37), 0·28 (−0·24, 0·79) and 0·56 (0·13, 0·96). Although the distributions are highly overlapping, the support at λ=1·51 does not contain the value of ρ=0. At the modal value λ=1·41, the residual skewness is practically zero, and the assumption of conditional normality is satisfied. Lower and higher values of λ generate slight negative and positive skewness, respectively (results not shown), and this results in the observed pattern of ρ.

The 95% confidence intervals for the MC estimates of posterior means indicate that the lengths of the chains (sample sizes) result in adequate accuracy in both analyses.

(d) Model comparison and posterior model checking

The Monte Carlo estimate of the pseudo marginal probabilities of the rabbit data under M and under M_λ are equal to −3930·7 and −3919·5, respectively. The corresponding figures for the pig data are −23 999·0 for M and −15 269·1 for M_λ. In both cases, the best fits correspond to the analysis on the transformed scale.

Figure 3 shows histograms of posterior realizations of the discrepancy measure (17) under M and under M_λ in the rabbit data and M and M_λ in the pig data, designed to test skewness at the level of the conditional distribution of the data. In both species, the distribution of (17) discloses negative skewness when the data are analysed in the original scale (first and third histograms from left), whereas under the transformed model, the mean residual skewness is zero (second and fourth histograms from left). This indicates that the transformation was successful in inducing normality at the level of the conditional distribution of the data.

Figure 3. Histograms of posterior predictive realization of , designed to test residual skewness for the rabbit data (first two from left) and pig data (last two from left). The first and third figures correspond to analyses on the original scale; the second and fourth to analyses on the transformed scale.

Since estimates of posterior means of σ_a*² and σ_p*² in the rabbit data were fairly small under M_λ, and sensitive to prior specifications, we investigated further whether heterogeneity could be accounted solely by systematic effects, without inclusion of random effects at the level of the environmental variance. In this model, labelled HM_λ, σ_a*²=0 and σ_p*²=0. The for HM_λ is −3927·3 indicating a poorer fit than the full model and thus support for σ_a*² and σ_p*². The corresponding figure for the pig data for HM_λ is −15 297·1, indicating also support for the full model.

Both data sets were also analysed using the classical infinitesimal model with homogeneous variance, conditional on the modal value of the Box–Cox parameter. The posterior mode of λ (95% posterior interval in brackets) is 1·41 (1·33,1·49), in the rabbit data and 1·26 (1·21,1·31) in the pig data. The under the classical infinitesimal model is −4009·5 in the rabbit data and −15 337·23 in the pig data. In the rabbit data, the result indicates that the global fit of the infinitesimal model is worse than that obtained fitting the genetically structured variance heterogeneity model M_λ (and also worse than model M). In the pig data, the global fit from the classical infinitesimal model is worse than that obtained fitting the genetically structured variance heterogeneity model M_λ, but better than the M model.

The overall conclusion from the analyses of both data sets is that the transformed model induces conditional normality of the sampling distribution of the data (as opposed to the untransformed model), and that the transformed model is best supported in terms of the overall fit relative to the untransformed model. Under the transformed model there is considerably weaker evidence in the rabbit data for a genetically structured variance heterogeneity model than that provided by the analysis in the original scale. In the pig data set, the transformed model provides evidence for the presence of genetic variation at the level of the environmental variance, but in contrast with the results of the analysis in the original scale, the correlation ρ between additive genetic values affecting mean and variance is strong and positive.

5. Discussion

This work has confirmed that inferences at the level of the residual variance are strongly affected by the scale on which the data are analysed. This is certainly not a surprising result, but it rather reflects the difficulties in learning about the state of nature via a statistical analysis of data using models that build on strong distributional assumptions. One option is to avoid introducing such assumptions and complex modelling, and instead use simple functions of the data specifically designed to uncover the structure associated with a genetic component at the level of variance. Uncertainty concerning the final result could be incorporated using resampling methods, such as permutation tests or the bootstrap. A limitation of this option is the difficulty of finding a measure that is free from artefacts and that signals unambiguously the existence of genetic variance on variance. Another limitation is whether the structure in the data related to the existence of genetic variance on variance is strong enough to be detected by these means. For example, we have attempted to regress sampling variance of repeated measurements within individuals in daughters on mothers and could not find a convincing pattern. A second option, followed in this work, is to use a fully parametric model that addresses specifically the putative presence of a genetic component affecting environmental variance. Support for the model is sought by comparing the quality of its fit with that of other models that neglect the presence of a genetic component on the environmental variance, using various criteria. An unambiguous favourable comparison is interpreted as evidence for the presence of a genetic component on the environmental variance. Spurious results can never be excluded using this approach either. The exercise reflects merely that, among the models considered, there is one that does better than the rest. It does not provide decisive evidence.

The results of the study differ in the rabbit and pig data. The transformed rabbit litter size data yields a Monte Carlo estimate of the posterior mean for the correlation (95% posterior interval in brackets) of 0·28 (−0·24, 0·79), compared to a value equal to −0·73 (−0·89, −0·50) obtained in the untransformed data. Further, the additive genetic variance at the level of the environmental variance based on the transformed data falls by a factor of 2·4, and its posterior distribution is considerably affected by prior assumptions. Thus, under the transformed scale, the support for a genetic component at the level of the environmental variance is markedly weaker. The picture is different for the pig litter size data. The correlation changes from −0·64 (−0·82, −0·45) to 0·70 (0·44,0·98) under the transformed model, the additive genetic variance at the level of the environmental variance falls by a factor of <2, and its posterior distribution is very little affected by prior input. On the basis of these results, the statistical evidence in favour of the genetically heterogeneous variance model is not disputed, but in the transformed model the correlation changes sign and its posterior distribution has positive support.

There has been some controversy concerning how inferences should be drawn using the Box–Cox family of transformations. Bickel & Doksum (1981) show that joint inferences of the transformation and the remaining parameters of a linear model lead to highly correlated estimates and the marginal variances of the latter are much larger compared to the conditional variances for fixed values of the transformation parameter. Hinkley & Runger (1984) and Box & Cox (1982) argue that inferences should be drawn conditional on a fixed value of the transformation parameter, since linear model parameters have meaning only with reference to a particular scale. This is an important point that argues against reporting inferences marginalizing over the distribution λ. The consequences on statements of posterior uncertainty fitting the model conditional on λ do not change the conclusions of our analysis. For example, for the rabbit litter size data, the 95% posterior interval of ρ changes from (−0·24, 0·79) for fixed λ, to (−0·36, 0·80) when the analysis is done marginalizing over λ. In the pig data, these figures are (0·44, 0·98) and (0·36, 0·99). In the context of prediction of future data, as in De Oliveira et al. (1997, 2002) and Lee et al. (2005), it may be sensible that posterior intervals of future data incorporate uncertainty about λ.

The problem of scale has clearly important implications in an animal breeding context and we consider two scenarios. First, assume that one is interested in selecting at the level of the mean only (via predicted additive genetic values) and that the ‘correct’ model is the transformed BCHET model. Here, the term ‘correct’ model is to be understood as the model under which the assumptions for the data in question hold. Does it make any difference whether predictions and subsequent selection decisions are obtained via the transformed BCHET or the untransformed HET model? For the pig data (the rabbit data show an almost identical pattern), the left-hand side of Figure 4 shows that there is a very close positive association between the posterior means of additive genetic values affecting the mean under both models, implying that approximately the same individuals are selected when ranked on the basis of either model. On the other hand, the right-hand side of Figure 4 discloses no association between additive genetic values affecting the variance under both models. Therefore, the efficiency of selection at the level of the variance will be impaired if selection operates on the basis of the ‘incorrect’ model.

Figure 4. Left: plot of the posterior means of the vector a under M (a_λ=1) versus the posterior means of the vector a under M_λ (a_λ=1·393). Right: plot of the posterior means of the vector a* under M (a_λ=1*) versus the posterior means of the vector a* under M_λ (a_λ=1·393*) in pig litter size data.

In the second scenario, assume again that the ‘correct’ model is BCHET and selection under this model (to change the mean, the variance, or both) leads to displacements at the level of additive genetic values affecting mean and variance. One may wish to know how changes in the transformed scale translate into changes in the original scale and report results in that metric. Under the BCHET model, genetic values operating on mean and variance are additive and the correlation between them has a clear interpretation. When transformed back to the original scale, genetic values are no longer additive and the correlation is difficult to define. However, the problem can be investigated by deriving the mean and variance of the conditional distribution of the data in the original scale, given knowledge of the transformed BCHET model. More specifically, consider the simplified version of the transformed BCHET model y^(λ)|λ, μ,μ*, a, a*~N (1μ+Za, exp (1μ*+Za*)). A Taylor expansion shows that the mean and variance in the original scale can be approximated by

(8)

and

(9)

respectively. For a given value of λ (and μ) one can study how the mean and variance in the original scale change due to displacements in (a,a*|ρ) on the transformed scale. Clearly, (a,a*|ρ) is non-linearly related to both (8) and (9) and Figure 5 illustrates some special cases. When λ=1·39, the variance in the original scale falls at a decreasing rate as a_i increases and a_i * is kept constant (implying ρ=0 in the transformed scale). The same scenario at a value of λ=0·8/0·0 shows that the variance in the original scale grows at a decreasing/increasing rate as a_i increases. On the other hand, selection for a_i * in the transformed scale, keeping a_i constant (implying again ρ=0 in the transformed scale), leads to a negative linear relationship between mean and variance in the original scale when λ=1·39. When λ=0·8 or λ=0·0, this relationship is linear and positive, respectively. When a_i and a_i * change simultaneously (results not shown), the relationship between mean and variance in the original scale is again highly dependent on λ and ρ. The overall conclusion is that it is not possible to make general statements about how changes taking place on the transformed scale translate into changes in the original metric.

In this work, the Box–Cox transformation was used to achieve normality (and not homogeneity of variance) at the level of the conditional distribution of the sampling model. The work can be extended in an obvious way by applying an unknown transformation to the additive genetic effects to induce normally at this level. For example under the homogeneous variance model, if data are conditionally normally distributed, but additive genetic effects a are not, then the marginal distribution of the data may display skewness. This can lead to spurious ρ when the (incorrect) HET model is fitted to the data. Other situations that can lead to misleading inferences include the presence of unaccounted mixtures, the wrong functional relationship between mean and variance, the wrong model for the sampling distribution of the data, or the presence of an unaccounted correlation between additive genetic and environmental effects. This makes it strikingly clear that an attempt to understand the state of nature via a statistical analysis of data must be regarded only as a first step until more fundamental knowledge becomes available.

Figure 5. Relationship between variance and mean (expressions (8) and (9)) when only a_i changes and a_i * is kept constant (top) and when a_i * changes and a_i is kept constant (bottom), for three values of λ.