Genetic heteroscedastic models for ordinal traits: application to sheep litter size.

BACKGROUND: Classical genetic canalization models, which accommodate the mean and variance of a trait separately, provide a flexible approach to take heteroscedasticity for continuous traits into account. However, this model is not appropriate for discrete traits. The aim of this work was to propose heteroscedastic threshold models suitable for the genetic analysis of ordinal data.
METHODS: In order to first fit the mean and variance of ordinal traits separately, we extended the classical threshold model (TM) for discrete data by introducing non-genetic and genetic factors of heterogeneity on the variance of its underlying variable, which leads to a homothetic threshold model HTM and its alternative parameterization HTM' in which the thresholds of different individuals are linked by a homothetic-translation. Relaxing the constraint between the thresholds led us to propose an independent threshold model ITM that was more flexible than HTM' but required the estimation of more parameters. TM, HTM and ITM models were applied to study 19,671 records on litter size in Romane sheep.
RESULTS: Both HTM and ITM were able to disentangle the link between the mean and variance that holds in the classical homoscedastic threshold model. The results obtained for the litter size of Romane ewes showed that the data was best fitted with HTM compared to ITM and TM. The correlations between the observed and predicted variances were equal to 0.6 and 0.2 for HTM and TM, respectively. These analyses showed the existence of a genetic component for the heterogeneity of litter size in sheep that was taken into account in HTM.
CONCLUSIONS: HTM is the most suitable model to study the variability of litter size in sheep. It accommodates both the mean and variance separately while requiring the estimation of only a few parameters.

Background

Robustness can be viewed as the ability to maintain a stable phenotype regardless of the environmental conditions. In line with this approach, Scheiner and Lyman [1] proposed a model which considers that the expression of a trait is controlled by two sets of genes. One set controls the level of the performance and the other set controls the environmental variability of the trait. Various authors developed heteroscedastic models (i.e. canalization models) to jointly estimate the influence of various factors on both the mean and the variability of a trait [2-5].

These models were implemented with frequentist [6, 7] or Bayesian methods [8, 9], and applied to a large variety of continuous traits (birth weight of piglets, birth weight of rabbits, weight of chickens, etc.).

The methods described above mainly target continuous traits. Theoretically, these linear models do not analyze discrete traits adequately because of their non-normal distributions. However, some previous empirical studies showed that the performance of linear and non-linear models was similar in terms of goodness-of-fit and predictive ability [10, 11] when analyzing the level of discrete livestock data. However, it is not known whether these properties are still true when dealing with the variability of discrete data.

The main objective of this research was to propose new models that take the heterogeneous individual variability of ordinal traits into account using derived threshold models. The models are described consecutively in the “Methods” section. “Heteroscedastic extensions of the threshold model” are presented in the first subsection, and a “Numerical application” of these models to litter size in Romane sheep is presented in the second subsection.

Methods

Heteroscedastic models for ordinal traits

Homothetic threshold models for ordinal traits

A simple way of studying environmental sensitivity of an ordinal trait is to combine the threshold model (TM) [12] (appropriate for the study of ordinal traits) and the canalization model (appropriate for continuous traits) [13]. Thus, heteroscedasticity is considered under the liability threshold model and we call this new model the homothetic threshold model or HTM:where is the normally distributed liability and is a set of thresholds, which transform this liability into the observed discrete variable in categories, is the overall mean, the vector of additive genetic effects that affect the mean of the liability, is the overall mean and is the vector of additive genetic effects that affect the variability of the liability. Several functions have been proposed in the literature (exponential, linear or square root functions [3]). The vectors of genetic values are assumed to be Gaussian:where denotes the Kronecker product, is the known additive genetic relationship matrix and the genetic matrix , where and are the variances of the additive genetic effects that affect the mean and the variability of the liability, and their coefficient of correlation.

A fully equivalent parameterization of HTM can be written. This parameterization or HTM’ model is reached by fixing the parameters of the liability and considering that the matrix of thresholds is affected by environmental and genetic effects i.e.:

Under this parameterization, the vector of thresholds is specific to each individual , with where is an overall set of thresholds, a subject-specific deviation from this overall set and a vector of size that reflects the specific deviation of the subject threshold from this overall set with:where

Under this parameterization, the correlation between and is set to 1 and the transformation of to (from individuals to is not a translation as in TM but a homothetic translation, which means that the ratio is constant for the whole population i.e.: HTM and HTM’ are equivalent models with the same number of parameters () being estimated. The relationship between and is:

The estimations obtained with HTM’ can be transformed into the estimations that would have been obtained with HTM. Overall, in the case of a linear function , the approach produces:whereas for the other parameter we have:

Finally, the genetic correlation is:

These formulas are detailed in “Appendix 1”. The heritability of the trait on the liability scale is estimated by the ratio of the animal variance that acts on the mean to the total variance of the underlying variable i.e.:

For instance in the case of a linear function :

Extension of the homothetic threshold model into the independent threshold model

A straightforward extension of the above HTM’ model can also be considered and consists in releasing the constraint of the homothetic-translation link between the thresholds of the individuals. The model therefore becomes an independent threshold model ITM, i.e.:where with being an overall set of thresholds, and a vector of size , which reflects the specific deviation of the subject threshold from this overall set with:

The total number of parameters in this model is . It should be noted that ITM cannot be parameterized like HTM using a fixed set of thresholds and a model on the liability . In this model, fixing the thresholds to modify the parameterization would imply that the underlying variable would no longer be Gaussian.

Due to the lack of constraints between the thresholds (except for the obvious relationship ), ITM fits a much larger range of situations than a homothetic threshold model, which assumes a strong constraint either on the homothetic translation relationship between thresholds (HTM’) or on the normality of the liability (HTM).

Interestingly, we found that another parameterization of ITM is a constrained multiple-trait model MTM (Appendix 2). Indeed, each observation of the ordinal trait can be considered as the result of exclusive binary traits:

The observed variable can thus be transformed into indicator functions which can be analyzed using a multivariate binary threshold model. Only indicator functions need to be considered in this MTM since the value of one of the indicator functions is fixed conditionally to the others. Each sub-model can be written as :

According to this model, the probability that an individual might achieve the observation or that the kth indicator function is equal to 1 is:and the probabilities should satisfy the relation This constraint must be added to the multivariate binary threshold model.

As for ITM, there are parameters to estimate i.e. thresholds and variance covariance parameters.

The relationship between the thresholds of MTM and ITM is as follows:

Under the latter models, one cannot really consider that the expression of the ordinal trait is controlled by two sets of genes. The combination of elementary genetic values for each individual results in the heteroscedasticity of the global trait. In this model, the heritability can be calculated for each class as:

The different parameterizations of the models are summarized in Table 1. HTM’ and ITM can be fitted using the ASReml software.

Table 1

Main characteristics of the different models

Models	Number of parameters	Thresholds	Underlying
		Population
Homoscedastic
TM	n	t^∗=t1∗,t2∗,t3∗,…,tn-1∗	Y∼N0,1
Heteroscedastic
HTM	n+2	t=0,1,t_3,…,t_n-1	Y∼Nμ+u,f^2(η,v)
HTM’	n+2	t^′=t1^′,t2^′,t3^′,…,tn-1^′	Y∼N0,1
ITM	n+2n-1/2	t^″=t1″,t2″,t3″,…,tn-1″	Y∼N0,1
		Animal i
Homoscedastic
TM	n	ti∗=t1i∗,t2i∗,t3i∗,…,tn-1i∗ t.i∗=t_i+μ_i+u_i	y_i=ε_i
Heteroscedastic
HTM	n+2	t_i=0,1t_3i,…,t_n-1i	y_i=μ_i+u_i+ε_i
HTM’	n+2	ti^′=t1i^′,t2i^′,t3i^′,…,tn-1i^′ t.i^′=ti^′+μi^′+ui^′+vi^′	y_i=ε_i
ITM	n+2n-1/2	ti″=t1i″,t2i″,t3i″,…,tn-1i″ t.i″=ti″+μi″+ui″	y_i=ε_i

Numerical application

Data

These models were applied to study litter size in Romane sheep. All data (litter size, environmental variation factors and pedigrees over ten generations) were extracted from the national database for genetic evaluation and research managed by the Institut de l’Elevage (French Livestock Institute) and the Centre de Traitement de l’Information Génétique (Genetic Information Processing Center, Jouy-en-Josas, France). Data consisted of records on litter size of 11,073 multiparous ewes from 23 flocks with at least four lambings from 2000 to 2013. The mean prolificacy (2.13) in this sample was high and consequently there was a high global variability in litter size. The main features of this dataset are in Table 2. To avoid extreme case problems [14], litter sizes of more than four lambs were pooled into a single class for all analyses, which resulted in four classes of litter size i.e. [1, 2, 3, 4+].

Table 2

Features of Romane sheep

Item
Number of records	19,671
Number of animals	11,073
Number of sires	1096
Number of dams	6550
Mean	2.13
Distribution of litter size (%)
Single	17.70
Twin	55.01
Triplet	24.18
Quadruplet and +	3.11

The fixed effects included in the analysis were selected to be as similar as possible to those found in others studies [15, 16]. This was done by comparing the likelihood-ratio of all the possible linear mixed models which included the following fixed effects: age of ewe (nine levels), time interval between lambing in months (seven levels), parity (seven levels), age of the ewe at first lambing in months (nine levels), season of lambing (two levels), as well as the random effects: flock and ewe effect. For this step of the analysis, the ewe effect was considered as random but the genetic relationships between individuals were ignored.

Estimation of parameters and comparison of methods

The three threshold models TM, HTM’ and ITM described above were applied to this data. The classical TM, used as a basis for comparison, included the previously selected fixed effects and the ewe as both a random genetic effect () and a random permanent environmental effect (i.e. ). The homothetic threshold model HTM’ considered that fixed, genetic and permanent effects were common for all thresholds and threshold-specific random genetic effects with constrained genetic (co)variance to ensure homothetic translation (i.e. ). Finally, an additional extension resulting in the independent threshold model ITM was fitted and it considered that the fixed and permanent effects were common for all thresholds and included a genetic effect that was specific to each threshold with no constrained genetic (co)variance (i.e. ).

ASReml software [17] was used to implement all the models.

Models were compared using the likelihood ratio test for nested models (TM vs. HTM’, LRT df = 2; TM vs. ITM, LRT df = 5) and the Bayesian Information Criterion (BIC) for non-nested models (, HTM’ vs. ITM, where is the restricted maximum likelihood of the model, the number of observations, and the number of fixed effects and covariance parameters) [18]. The models were also compared based on their ability to predict new records as follows: 75 % of records were used to estimate parameters that were used to predict the remaining 25 %. Five replicates of this design were sampled at random. To evaluate the predictive ability of the models, MSEP (mean square error of prediction) statistic was computed as:where and correspond to the observed and predicted litter size, is the number of data in the validation subset and is the subset number.

For ewes with observations, we compared the raw means and variances of their litter sizes to the corresponding predicted values on the observable scale obtained with the different models. The raw mean and variance were calculated for each individual from the observed data. The predicted values were calculated as:and where was the probability that the litter size of a ewe was . These probabilities were calculated using with as above and ( corresponds to or depending on the model). To measure the ability of the models to break the link between and we used the following criterion:which ranges from 0 to 1. The criterion was 0 if which means that a single value of the mean has a single value of the variance. The criterion was 1 if which means that a single value of the mean has different variance values. Because the density associated to the variance is unknown, data were ordered in ascending order of litter size and were discretized into subgroups, and the variance of y in each subgroup was calculated in order to calculate:

The subgroup number was defined in order to get a sufficient amount of data in each subgroup to obtain a good estimation of the subgroup variance. The intensity of the dispersion of each model was also assessed by calculating the correlation, either for the mean or the variance, between the observed values and the predicted values on the observable scale.

Results

Comparison of models

The values of the criteria (LRT, BIC, MSE and MSEP) used to compare the adjustment quality of the models are in Tables 3 and 4.

Table 3

Values obtained from likelihood-ratio, BIC tests, B criteria, and mean squared error (MSE) between observed and fitted values in the different models

Models/tests	Likelihood-ratio (LRT)	BIC test	B	MSE of mean	MSE of variance
TM	TM versus HTM': ≫11(Reject H0 P value <0.001)TM versus ITM: ≫16(Reject H0 P value <0.001)	54,005	0.013	0.07	0.13
HTM’		49,517	0.51	0.07	0.10
ITM		50,663	0.42	0.07	0.10

TM threshold model, HTM’ homothetic threshold model, ITM independent threshold model, B B criteria

Table 4

Values obtained from mean squared error of prediction (MSEP) and B criteria on the validation subset in the different models

Models/tests	MSEP	B criteria
Threshold model (TM)	0.473	0.005
Homothetic threshold model (HTM’)	0.471	0.363
Independent threshold model (ITM)	0.481	0.371

The results showed that the homoscedastic model should be rejected and that the models that account for genetic heterogeneity of litter size distribution should be favored. Among the heteroscedastic models, BIC values were lower with HTM’ than with ITM. MSE estimates of the mean were similar for the three models, whereas MSE estimates of the variance differed and were higher with TM than with HTM’ (0.13 with TM vs. 0.10 with HTM’). The differences between MSEP values are small, which is explained by the fact that this parameter can estimate prediction error for the level of the trait only and not for its variability. To evaluate the ability of the models to predict variance, we calculated the B criteria on the validation subset by discretizing the whole data into 200 adjacent subgroups. The results are in Table 4. Based on the B criteria, HTM’ and ITM were more efficient than TM (0.36 for HTM’, 0.37 for ITM vs. 0.005 for TM).

Homoscedastic threshold model used as reference (TM)

With TM, the genetic additive variance of the liability was equal to and the variance of the permanent environmental effect was equal to , which resulted in a heritability of 0.07. The correlation between the predicted and observed means was very high () (Table 5), but it was low between the predicted and observed variances (). In spite of the curvilinearity between the predicted means and variances, the correlation between these variables was high () and close to 1, and the parameter was equal to 0.013 (Table 3), which reveals the one-to-one relationship between these variables when this homoscedastic threshold model was used as shown on the scatterplot of Fig. 1a.

Table 5

Correlations between the predicted and observed means as well as the predicted observed variance with the different models

Model	Predicted mean versus observed mean	Predicted variance versus observed variance
TM	0.80	0.21
HTM’	0.79	0.61
ITM	0.78	0.59

TM threshold model, HTM’ homothetic threshold model, ITM independent threshold model

Fig. 1

Scatterplots of estimated mean vs. estimated variance across the observed range for the three models. a TM model. b HTM' model. c ITM model

Homothetic threshold model (HTM’)

The genetic additive variances obtained with HTM’ are in Table 6. By design, the correlation between the genetic values associated to thresholds () was equal to 1 in this model and the variance was equal to 0. Moreover, and as expected due to the homothetic relationship, we observed an increase of compared to . The variances were equal to 0.04, 0.05 and 0.20 for , respectively. The variance of the permanent environmental effect was equal to 0.015.

Table 6

Genetic variance and correlation between thresholds for Romane sheep using the HTM’ model

Thresholds/genetic variance	σu^′.12	σν^′.22	σν^′.32
	0.04	0.13	0.13
		0.05	1
			0.20

By transforming the results of the HTM’ parameterization to those of HTM, we could estimate the variance of the genetic effects that affect the mean of the liability as well as the correlation between these effects. The values found for the mean of the liability, the scaling function and the correlation between effects were , and , respectively. The variance of the permanent effect was whereas, by design, the mean of the residual variance was . Since the model fitted for the environmental variance did not include fixed effects but only a general mean and a genetic effect, the heritability of the liability mean could be estimated independently of a given environment and was as follows:

The correlation between observed and predicted means was high () and not different from the value obtained with TM. This contrasts with the correlation between predicted and observed variances (). Although the prediction of the variances was still lower than that of the means, it was much higher than with TM. HTM’ could therefore break the link between the predicted means and variances and the B parameter was equal to 0.51 (Table 3) even if the correlation between these variables was still high (). The scatterplot in Fig. 1b shows this variability of the predicted variance for a given predicted mean.

In this homothetic translation model, the gap between the two first thresholds is an indirect measure of its variability for each individual. As an example, we estimated the distribution of the litter size of two individuals having the same average litter size but different variances (Table 7). For these individuals, the values of the mean and variance estimated with HTM’ are close to the observed values of the mean and variance. The ratio is constant for the two individuals, and the difference in percentage of grand-offspring born as quadruplets or more can be large (0 vs. 6 %).

Table 7

Observed values (prolificacy, variance) and estimated values with HTM’ [thresholds, different gap between the first thresholds (), ratio (R), prolificacy, variance, and percentage of the simple (% LS 1), double (% LS 2), triplet (% LS 3) and quadruplet (% LS 4)] for two individuals

Observed values			Estimated values with HTM’ model
Ind	Prolificacy	Variance	τ1^′	τ2^′	τ3^′	Δ_1-2	R	Prolificacy	Variance	% LS1	% LS2	% LS3	% LS4+
1	2.00	0.50	−0.63	0.92	2.16	1.55	1.8	2.09	0.54	16	53	25	6
2	2.00	0.33	−0.89	0.86	2.28	1.75	1.8	2.02	0.38	21	64	15	0

Independent threshold model (ITM)

The genetic variances and correlations obtained with ITM are in Table 8. Genetic variances were small, but larger for the second and third thresholds than for the first threshold. Genetic correlations were moderate () and it was highest between and . Using the transformation formulas described in “Appendix 2”, it was possible to estimate the genetic variance associated with twin lambing (), the corresponding heritability () as well as the genetic correlation between the values for the first and the second classes of litter size (). The correlations between the predicted and observed means () and between the predicted and observed variances () were close to those obtained with HTM (Table 8). The parameter (, Table 3), the correlation between the predicted means and variances () and the scatter plot (Fig. 1c) were also similar to those obtained with HTM’. The Spearman rank correlations between the thresholds of HTM’ and ITM were high (Table 9).

Table 8

Genetic variance and correlation between thresholds for Romane sheep using the ITM model

Thresholds/genetic variance	σν^″.12	σν^″.22	σν^″.32
	0.035	0.631	0.430
		0.138	0.921
			0.176

Table 9

Spearman rank correlations between the HTM’ and ITM models

First threshold	0.988
Second threshold	0.989
Third threshold	0.988

Discussion

A genetic heteroscedastic model was previously proposed [6] to analyze categorical variables and break down the link between mean and variance that holds with the classical homoscedastic threshold model. This model originated from the combination of two base models: the homoscedastic threshold model for categorical data [19] and the structural heteroscedastic model for continuous variables [2]. In this paper, we studied further extensions of this model by analyzing the relationship that exists between categories and in particular by releasing the main constraint.

In this paper, we propose two kinds of heteroscedastic models, both with different parameterizations. The first type (HTM/HTM’) was fully derived from the canalization model proposed for continuous traits [2] with which it was easy to draw comparisons. The second type consisted in relaxing the relationship that exists between thresholds (translation in TM or homothetic translation in HTM) and led to a model that considers each threshold independently i.e. ITM. We demonstrated that this model was equivalent to a constrained multiple-trait model MTM on each class of the ordinal trait. This second type of model is much more flexible than HTM and its application in selection would be of particular interest because it allows for each class to be controlled through a multiple selection index [20] with specific environmental effects and economic weights for each class. Modeling each class of the trait to control its distribution or modeling only the most important class are trivial issues and have already been considered [21], however the correct multiple trait model must include a constraint on the sum of the probabilities of the binary variables. In spite of their potential interest, the use of this kind of model (ITM or MTM) for genetic evaluation and selection would be limited since they require a larger number of parameters to be estimated than HTM or HTM’ and would probably show convergence problems when the number of classes exceeds three. Furthermore, their biological interpretation is difficult; for instance we could not link the parameters of ITM or MTM to the biological model of Schneider and Lyman [22], which assumes the existence of different genes that control the mean and the variance of a trait.

We compared TM, HTM’ and ITM through the analysis of data on litter size from both Romane and Rouge de l’Ouest sheep (the latter results are not shown). The low heritability of the liability to litter size obtained with an homoscedastic threshold model (i.e. TM) agreed with previous studies that used TM [23-25]. However, LRT, BIC, and MSE parameters showed that this homoscedastic model was less successful than both heteroscedastic models. Consequently, the distribution of sheep litter size shows heteroscedasticity of genetic origin as revealed by both the homothetic HTM’ and independent threshold ITM models.

Taking the genetic heteroscedasticity in HTM’ and the subsequent transformation of the parameters to those of HTM into account did not substantially change the heritability of the liability to litter size compared to TM. However, with HTM, the additive genetic variance of the variability of the trait was smaller than the corresponding variance of its mean. This is in line with results in the literature that showed that when additive genetic variances of the mean and the variability of a trait were jointly estimated, their values were smaller for the variability than for the mean [6, 9, 26–31]. This seems to be a common rule regardless of the methodology used for such estimations, the trait or the species. In contrast, the genetic correlation between the mean of a trait and its variability is highly variable and covers the whole range of its variation interval. For instance, in a mouse population, genetic correlations of −0.93 and 0.97 were reported for litter size and individual birth weight, respectively [29]. In our study, we found a negative genetic correlation between the mean and the variability of sheep litter size with HTM’ (). As highlighted by Yang et al. [32], this parameter might be affected by an artefact due to the scale of measurement or skewness of the data, although with the threshold model no such skewness affected this correlation. However, a study on another French sheep breed (Rouge de l’Ouest) for which the mean prolificacy (prolificacy = 1.80; LS3+ = 12.5 %) is lower than in Romane sheep (results not shown) reported a similar value i.e. −0.25.

Relaxing the homothetic constraint between the thresholds of HTM’ did not improve the fit of the model or the predictive ability of the data as might be expected, but instead increased the BIC values due to the larger number of parameters to estimate in ITM. Moreover, although the constraint on the genetic correlation between thresholds was relaxed in ITM, the correlation between the second and third thresholds remained close to 1 and the correlation between the first and second thresholds was high; in addition, the Spearman rank correlations between the thresholds of HTM’ and ITM were close to 1 (Table 9). This suggests that the true nature of the link between the thresholds is almost a homothetic translation in sheep and that the distributions of litter size classes are not independent but linked by a strong common law. It also means that the concept of sensitivity thresholds, which transform individual Gaussian variables that differ in mean and variance into the observed categories, might be relevant.

Nonetheless in our study, only genetic effects were considered as threshold-specific in ITM in order to compare the different models using LRT or BIC criteria. Of course, threshold-specific fixed effects and permanent effects could be included in ITM. In the same way, threshold-specific fixed effects and permanent effects could be included in HTM’. Unfortunately, HTM with a threshold-specific permanent effect did not converge and threshold-specific fixed factors with proportional effects cannot be implemented in ASReml. It may be possible to estimate fixed effects in the residual variance using an iterative process that tests the influence of each effect on the variance of the model. We did not test such algorithms on our data and it would probably be a slow and difficult process. Another solution would be to modify existing software for canalizing selection of a continuous trait in order to take into account the discrete characteristics of the data.

Using ITM, Martin et al. [33] showed that the polymorphism of a single gene affected not only the mean prolificacy in Lacaune sheep, but also modified its variability by changing the gap between the thresholds in a homothetic way. A superiority of the homothetic model (HTM’), as well as a similar genetic variability of the liability to litter size were found for Rouge de l’Ouest sheep (results not shown), thus covering the whole range of prolificacy in breeds for which canalization might be of interest. By design, the gap between the first two thresholds estimated for each individual in HTM’ is inversely proportional to the individual residual variance of HTM (Fig. 2). This criterion has been used to compare the variability of about 50 French sheep breeds [34].

Fig. 2

Scatterplot of estimates of the gap between the threshold and the mean predicted by the HTM model

Conclusions

This study shows that we can model the variabilities (mean and variance) of a discrete trait using a suitable model, namely the HTM’ model. Thus, the canalization of discrete traits is possible and the gap between the first two thresholds estimated for each individual in this model could be used to estimate the breeding value for the variability of litter size independently from that of the mean.