J. Dairy Sci. 87:699-705
© American Dairy Science Association, 2004.
Random Regression Test-Day Models with Residuals Following a Student’s-t Distribution
J. Jamrozik1,
I. Strandén2 and
L. R. Schaeffer1
1 Centre for Genetic Improvement of Livestock, Department of Animal and Poultry Science, University of Guelph, Guelph, ON, N1G 2W1, Canada
2 MTT Agrifood Research Finland, Animal Production Research FIN-31600 Jokioinen, Finland
Corresponding author: J. Jamrozik; e-mail: jamrozik{at}sherlock.aps.uoguelph.ca.
 |
ABSTRACT
|
|---|
First-lactation milk yield test-day records of Canadian Holsteins were analyzed by single-trait random regression test-day models that assumed normal or Student’s-t distribution for residuals. Objectives were to test the performance of the robust statistical models that use heavy-tailed distributions for the residual effect. Models fitted were: Gaussian, Student’s-t, and Student’s-t with fixed number of degrees of freedom (equal to 5, 15, 30, 100 or 1000) for the t distribution. Bayesian methods with Gibbs sampling were used to make inferences about overall model plausibility through Bayes factors, posterior means for covariance components, estimated breeding values for regression coefficients, solutions for permanent environmental regressions, and residuals of the models. Bayes factors favored Student’s-t model with the posterior mean of degrees of freedom equal to 2.4 over all other models, indicating very strong departure from normality. Number of outliers in Student’s-t model was reduced by 35% in comparison with the Gaussian model. Differences in covariance components for regression coefficients between models were small, and rankings of animals based on additive genetic merit for the first two regression coefficients (total yield and persistency) were similar. Results from the Gaussian and Student’s-t models with fixed degrees of freedom become more alike (smaller departures from normality for Student’s-t models) with increasing number of degrees of freedom for the t-distributions. For any pair of Student’s-t models, the one with the smaller number of degrees of freedom for the t-distribution was shown to be superior. Similarly, number of outliers increased with increasing degrees of freedom for the t distribution.
Key Words: random regression model • test-day data • Student’s-t distribution
Abbreviation key: ASD = age at calving by season of calving by DIM interval, BF = Bayes factor, GS = Gibbs sampling, PE = permanent environment, RR = random regression, TD = test day
|
INTRODUCTION
|
|---|
Outliers or abnormal test-day (TD) yields are common in dairy cattle. Reasons for atypical, infrequent observations include measurement errors, sickness or preferential treatment of cows on a given TD, short-term changes in herd environment, and mismanagement of the data. Wiggans et al. (2003) reported 1.8% abnormal milk yields (defined as <60 or >150% of predicted TD yield) in over 93 million TD records of Holsteins in the United States. Such abnormal measurements may adversely impact genetic evaluation if not accounted for properly. Human errors or manipulation (for example preferential treatment of dairy cattle) may unduly increase genetic merit of treated animals and their relatives leading to incorrect breeding decisions.
Several methods to detect and correct outliers prior to data analysis have been proposed. Other approaches use robust techniques to accommodate drastic errors contaminating the data in statistical models. Robust methods include modeling random effects (the residual term in particular) with heavy tailed distributions, like the Student’s-t (Lange et al., 1989). The t-distribution has thicker tails than the normal distribution and more variation (including abnormal yields) in the data is allowed. Influence of outliers may therefore be less profound and inferences about parameters more correct.
Random regression (RR) TD models have been based on normality assumptions for the random effects (Jamrozik and Schaeffer, 1997). Models with Student’s-t distributions for residuals have been developed for analysis of 305-d yields (Strandén and Gianola, 1997, 1999) and were compared with the Gaussian model for several scenarios. Results indicated that Student’s-t models outperformed their Gaussian counterparts in terms of model comparison statistics with Bayes factors, and analysis of outliers (observation that gave extreme residual terms in the model).
Objectives of this research were to compare Gaussian random regression TD models with several models that used Student’s-t distribution for residual terms. Comparisons included Bayes factors, estimates of covariance components and breeding values, solutions for permanent environmental (PE) regression coefficients, and distribution of outliers. Analyses were conducted via Bayesian methods employing Gibbs sampling (GS).
|
MATERIALS AND METHODS
|
|---|
Data
Data were first-lactation TD milk yields of Canadian Holsteins used previously in the international TD model study for 4 countries: Australia, Canada, Italy, and New Zealand (Jamrozik et al., 2002). The TD data comprised 8,839,030 records on 1,075,516 cows in 14,486 herds. Cows calved between January 1, 1990, and December 31, 1997. Overall edits were: DIM between 5 and 305, age at calving between 18 and 38 mo, and daily milk yield between 0.1 and 100 kg. Only records corresponding to 2 times a day milking were used. A small subset of the data including 120,446 records on 14,798 cows in 150 herds was created. The following steps were applied to select the data:
Fifty sires with daughters in Canada and in at least one other country (Australia, Italy, and New Zealand) were chosen in random; 150 herds with at least one daughter of selected sires were chosen in random; TD records on all cows in selected herds were kept.
Total number of animals in the analyses was 25,087. Contemporary groups were herd-TD classes with 6107 levels. Two seasons of calving (March–August, September–February) and 3 age at calving classes (<26, 26 to 30, >30 mo) were used to define 6 age x season groups. Days in milk on TD were divided into 29 classes: 5 to 20, 291 to 305, and 10-d intervals from DIM 21 to 290. Age at calving by season of calving by DIM interval (ASD) classification included 174 levels.
Models
The overall single-trait random regression model was a 2-stage linear hierarchical model (Jamrozik et al., 2001) that can be written as:
where yi is a vector of ni TD records of cow i (i=1,2,...,N) taken at known times ti, h is a vector of herd-test date effects, b is a vector of ASD effects, matrices Hi and Xi relate cow records to appropriate elements in h and b, respectively. f (wi, ti) is an expected trajectory (lactation curve) for ith cow, with cow-specific parameters wi that was modeled using Legendre polynomials of order 4, 
i is a vector of ni residuals associated with yi. In the second stage model, ui is a vector of additive genetic effects of cow i on wi, and ei represents a vector of second stage residuals (PE effect).
Student’s-t model.
The first stage conditional distribution for Student’s-t model was
where tni is ni dimensional t-distribution with v degrees of freedom. The second stage conditional distribution was wi | ui, E ~ N [ui, E].
Prior distributions for the parameters were:
where A is an additive genetic relationship matrix between individuals, and G is the additive genetic covariance matrix between elements of ui,
where 
and s2 are parameters of inverted chi-square distribution, g (e) and G0 (E0) are hyper-parameters of the inverted Wishart distributions. Prior distribution for degrees of freedom parameter () was a uniform distribution, continuous in the interval [2.0, 30.0].
Gaussian model.
The first-stage conditional distribution for the Gaussian model was
All other distributional assumptions (excluding degrees of freedom) were the same as for the Student’s-t model.
Seven different models were fitted and they were:
Gaussian — model with normal distribution of residuals at the first stage,
t — model with Student’s-t distribution of residuals at the first stage, and
t(i) (i = 5, 15, 30, 100, 1000) — models with Student’s-t distribution of residuals at the first stage, where i = the degrees of freedom, which was assumed to be known.
Model Fitting
Analyses were conducted via Gibbs sampling with Metropolis-Hastings step for degrees of freedom () in the t model. Detailed sampling algorithms for the Gaussian model can be found in Jamrozik et al. (2001). Strandén and Gianola (1999) gave corresponding algorithms for Student’s-t models. Prior values for covariance components with small degrees of freedom (24, 24, and 50 for genetic, PE, and residual covariances, respectively) were the same for all models and were taken from Jamrozik et al. (2002). Starting values for dispersion parameters were the same as corresponding priors and starting value for degrees of freedom of Student’s-t distribution was equal to 10. Proposal distribution for degrees of freedom of the t-distribution was N[vcurrent, vcurrent/10 ], where vcurrent is the current value of v. Samples (520,000 in total) were generated for each model with 20,000 iterations as a burn-in. Convergence of GS was monitored by visual inspection of selected parameters. Posterior means and standard deviations were calculated for all parameters in each model.
Model Comparisons
Bayes factors.
Models were compared using Bayes factors (BF) (Kass and Raftery, 1995). The BF to contrast model M0 (with parameters 
0) against model M1 (with parameters 1) is defined as B01 = p(y | M = M0) / p(y | M = M1), where p(y | M = Mk) = 
p(y | k, M = Mk) p(k | M = Mk) dk denotes an integrated (marginal) likelihood, and p(k | M = Mk) is the prior density under model k. If the prior probabilities of each model being true are equal, then the BF gives the ratio between the posterior probabilities of each pair of models. Values of B01 greater than 150 (or log of B01 is greater than 10) indicate very strong evidence in favor of model M0. The marginal density p(y | M = Mk) was estimated by the harmonic mean of the likelihood values from the Gibbs chain as 
, where 
(i =1, 2, ..., m) were draws obtained from the posterior distribution under model k.
Outliers.
A TD observation was defined as an outlier in a given model when the absolute value of the posterior mean of corresponding residual effect divided by the posterior standard deviation of this residual was larger than 3.
|
RESULTS
|
|---|
Numbers of independent samples (Geyer, 1992) in the t model ranged from 815 to 2048 and from 1184 to 2617 for genetic and PE covariance components, respectively. There were 40,620 independent samples for residual variance and 92,022 independent samples for degrees of freedom for the t model. The Gaussian model gave similar numbers of independent samples.
The acceptance rate for the Metropolis-Hastings step in GS for the t model was 24%. The posterior mean degree of freedom for the t model was 2.41, with standard deviation of 0.04. The values of these estimates were invariant (results not shown) to the starting values for degrees of freedom used for GS scheme. Figure 1
shows posterior distribution of degrees of freedom for the t model. The distribution was extremely sharp, as indicated by posterior standard deviation. All values for degrees of freedom were concentrated in the interval from 2.20 to 2.60.
View larger version (13K):
[in this window]
[in a new window]
|
Figure 1. Estimated posterior marginal distribution of degrees of freedom (df) for Student’s-t (t) model.
|
|
Estimates of marginal likelihoods for Gaussian, t, and t(i) models (for i = 5, 15,30,100,1000) are shown in Table 1. Bayes factors (natural logarithm of the Bayes factor is the difference between logarithms of appropriate marginal likelihoods) indicated very strong superiority of the t model over the Gaussian model and all models with constant degrees of freedom for the t-distribution. Differences between the Gaussian model and the Student’s-t models decreased with increased number of degrees of freedom for the t-distribution. Increasing degrees of freedom in t(i) models resulted in decreased overall plausibility of the model. Student’s-t model with 1000 degrees of freedom, however, still differed very strongly from the Gaussian model. Estimates of posterior means for residual variance (Table 1) showed very similar pattern in terms of model comparisons. Residual variance of the t model was 35% of the Gaussian’s residual variance. Estimates of residual variance for t(i) models increased with i, values for t(1000) and the Gausian models were practically identical.
View this table:
[in this window]
[in a new window]
|
Table 1. Posterior means (SD in brackets) of residual variance ( 2) and estimates of marginal likelihood (log p(y | M = Mi)), by model.
|
|
Tables 2 and 3 give posterior means of covariance components for genetic and PE regression coefficient, respectively, for Gaussian and the t models. Estimates from both models showed a high degree of similarity with slightly bigger differences for PE components. Heritabilities of the intercept of the Legendre curve were 0.46 and 0.48 for the Gaussian and t models, respectively. Corresponding values for the linear term were 0.48 and 0.51. Correlations between the first two PE coefficients were the same for both models, genetic correlation decreased (from -0.11 to -0.19) when the t model was fitted. Plots of daily realizations of genetic and PE variances for both models are given in Figure 2.
View this table:
[in this window]
[in a new window]
|
Table 2. Posterior means (SD in brackets) of covariance components for genetic regression coefficients (u0, u1, u2, u3, and u4), for Gaussian and Student’s-t (t) models.
|
|
View this table:
[in this window]
[in a new window]
|
Table 3. Posterior means (SD in brackets) of covariance components for permanent environmental (PE) regression coefficients (e0, e1, e2, e3, e4), for Gaussian and Student’s-t (t) models.
|
|
View larger version (17K):
[in this window]
[in a new window]
|
Figure 2. Daily genetic and permanent environmental (PE) variances for Gaussian (normal) and Student’s-t (t) models.
|
|
Estimates of linear regressions and correlations of posterior means for regression coefficients of Gaussian on Student’s-t models are in Tables 4 and 5 for genetic (25,087 animals) and PE (14,798 cows) components, respectively. Similarity of means for regression coefficient between Gaussian and various Student’s-t models was good, in general, with slightly better agreement for additive genetic effects. Lower order coefficients showed bigger differences, especially for smaller number of degrees of freedom for t distributions. Models t(100) and t(1000) gave virtually the same posterior means for regression coefficients as the Gaussian model.
View this table:
[in this window]
[in a new window]
|
Table 4. Comparison of estimated breeding values for genetic regression coefficients (u0, u1, u2, u3, and u4), between Gaussian and Student-t models.
|
|
View this table:
[in this window]
[in a new window]
|
Table 5. Comparison of solutions for permanent environmental (PE) regression coefficients (e0, e1, e2, e3, e4), between Gaussian and Student’s-t models.
|
|
Figure 3 shows a plot of posterior means of ASD effect for Gaussian and t models for two extreme age of calving—season of calving categories (season March–August for the youngest vs. season September–February for the oldest cows). The shape of curves for competing models was the same for both classes. Curves pertaining to the t model, however, were higher in overall level compared with their corresponding Gaussian counterparts. Similar pattern was observed for the remaining ASD classes (results not shown).
View larger version (19K):
[in this window]
[in a new window]
|
Figure 3. Age at calving (A) by season of calving (S) effects for Gaussian (normal) and Student’s-t (t) models. A1 = <26 months, A3 = >40 months, S1 = March–August, S2 = September–February.
|
|
Numbers of outliers (overall, with plus sign and with minus sign) for all analyzed models are presented in Table 6. The Gaussian model exhibited 3.7% of the total number of records with large relative values of residuals. Model t gave approximately 35% less observations that showed outlier properties. As number of degrees of freedom in t model increased so did the number of detected outliers. Models t(i) for i >30 and the Gaussian model gave roughly the same number of outliers. Distributions of outliers for all models were skewed, with more observations having negative values of residuals. Distributions of outliers by DIM intervals for Gaussian and t models (Figure 4) showed similar properties, with the majority of observations in the middle of lactation and relatively small frequency of outliers at the edges of the trajectory. Model t gave a smaller number of outliers in every DIM class than the Gaussian model.
|
DISCUSSION
|
|---|
Models for 305-d yields with Student’s-t distribution for residuals have been shown to outperform models with normality assumption for residual effect (Strandén and Gianola, 1997, 1999). Random regression TD models showed similar behavior, as indicated by the results of this study. The estimate of posterior mean for number of degrees of freedom was extremely low. This, however, was close to the finding of Strandén and Gianola (1997), in which the corresponding estimates ranged from 5.30 to 5.96. Test-day observations might be more prone to aberration, thus thicker tailed distribution for residuals may be more useful to accommodate such data. Attempts were undertaken to fit the Student’s-t model with the uniform discrete prior distribution for the degrees of freedom. Results (not shown) indicated very good agreement with the ‘continuous prior’ model as presented in this paper. Posterior mean of degrees of freedom was in particular also very small.
Student’s-t distribution for residuals is equivalent to the mixture of normal distributions with the ‘scaled’ residual variance equal to 
, where 
for the ith cow. Estimates of residual variances for the Student’s-t models were significantly smaller than the Gaussian residual variance (at least for lower number of degrees of freedom). ‘Scaling’ factors appropriately adjusted residual variance in accordance with cow specific 
. Separate (conditionally independent) t-distributions were applied to different cows in our model. The cluster t-model allowed for different scaling parameters for records on different cows.
Changes in the value of residual variance between Gaussian and the t model (Table 1
) resulted in proportional (on average) changes in solutions for fixed effects. This scaling effect could (at least partially) explain systematic differences between estimates for ASD effects between models (Figure 3).
Animal genetic and PE effect were modeled in a traditional way using normality assumptions. Further studies could assume heavy-tailed distributions for these effects also. However, it may be necessary to study the existence of heterogeneous variance in the data because a heavy-tailed distribution for the genetic effect is similar to a heterogeneous genetic variance model.
Bayes factors favored Student’s-t models over the Gaussian model. Similarly, for any pair of Student’s-t models, the one with the smaller number of degrees of freedom for the t distribution was shown to be superior. Differences between Gaussian and the t model were very big in terms of BF. Causes might include preferential treatment of some animals (daughters of sires used internationally).
Student’s-t distribution with >30 degrees of freedom usually has been treated as the normal distribution. This approach has not been confirmed by the values of Bayes factors estimated in this study. Models with the Gaussian and Student’s-t distribution with 1000 degrees of freedom were not the same in terms of Bayes factors. The relatively large number of observations (N) could play an important role. As N goes to infinity, the posterior probability of selecting the correct model by Bayes factor goes to 1 (O’Hagan, 1994). It should be noted also that the model with fixed number of degrees of freedom (models t(i)) and the t model are conceptually not the same.
Superiority of heavy-tailed models over the Gaussian one was confirmed by the outlier analysis. The proportion of outlier observations in the RR TD model with normality assumptions was larger than reported by Strandén and Gianola (1997) for the 305-d model (3.7 vs. 1.4%) even for slightly more restrictive definition of outliers in the current study (3 vs. 2 posterior standard deviations). Usefulness of Student’s-t distribution for analysis of TD data with RR model could potentially be greater than in the case of the model for total yields in lactation. Reduction in the number of outliers using the t-distribution was slightly larger in Strandén and Gianola (1997) compared with the current study.
Estimated breeding values for the 305-d yield in lactation gave correlation of 0.991 between the Gaussian and the t models. Rankings of animals would practically be the same by both models for this trait. Persistency of lactation (defined as the linear term coefficient) showed a lower degree of correlation between these 2 models. In practical terms, significant reranking of animals between Gaussian and t-distributions for this trait could be expected. Bigger differences between models were detected for estimates of regression coefficients for PE effect. Therefore, reducing the influence of outliers through robust distribution for residuals had a more profound effect on environmental parameters.
Outliers increase the variability of records within the contemporary group. Therefore, several models for analysis of dairy production traits use heterogeneous variance adjustment to control (at least partially) the effect of outliers. Within herd-test day correction factors also are used in a routine genetic evaluation model (Canadian TD model) in Canada (Schaeffer et al., 2000). All records in a contemporary group with extreme residual variance are subject to adjustment in this approach and not just the outliers. Student’s-t model, on the other hand, scales records on a cow basis. No precorrection for heterogeneous variances was done in the current study. Differences between Student’s-t models and the Gaussian model could have been possibly smaller if such preadjustment had been done.
|
CONCLUSIONS
|
|---|
Superiority of modeling residual effect in RR TD model with Student’s-t distribution over the usually assumed normality scenario was documented through Bayes factors and analysis of outliers. The best model was the Student’s-t model with the posterior mean of degrees of freedom equal to 2.4. This indicated very strong departure from normality for the residuals of TD data. Estimated breeding values for economically important traits derived from regression coefficients were affected to a lesser extent by using heavy tailed distributions for residuals. Fitting the model with Student t-distribution required Bayesian tools with GS. For large populations of dairy animals this seems to be beyond the current technical possibilities.
|
ACKNOWLEDGEMENTS
|
|---|
Data were provided by the Canadian Dairy Network of Guelph, Ontario, Canada. The authors are grateful to the Ontario Ministry of Agriculture, Food and Rural Affairs, DairyGen, and the Natural Science and Engineering Research Council for their financial support.
Received for publication June 19, 2003.
Accepted for publication August 10, 2003.
|
REFERENCES
|
|---|
Geyer, C. J. 1992. Practical Markov Chain Monte Carlo. Statist. Sci. 7:473–483.
Jamrozik, J., D. Gianola, and L. R. Schaeffer. 2001. Bayesian estimation of genetic parameters for test day records in dairy cattle using linear hierarchical models. Livest. Prod. Sci. 71:223–240.
Jamrozik, J., and L. R. Schaeffer. 1997. Estimates of genetic parameters for a test day model with random regressions for yield traits of first lactation Holsteins. J. Dairy Sci. 80:762–770.[Abstract/Free Full Text]
Jamrozik, J., L. R. Schaeffer, and K. A. Weigel. 2002. Estimates of genetic parameters for single- and multiple-country test-day models. J. Dairy Sci. 85:3131–3141.[Abstract/Free Full Text]
Kass, R. E., and A. E. Raftery. 1995. Bayes factors. J. Am. Stat. Assoc. 90:773–795.
Lange, K. L., R. J. A. Little, and J. M. G. Taylor. 1989. Robust statistical modeling using the t-distribution. J. Am. Stat. Assoc. 84:881–896.
O’Hagan, A. 1994. Kendall’s Advanced Theory of Statistics Vol 2B: Bayesian Inference. Oxford University Press, New York, NY.
Schaeffer, L. R., J. Jamrozik, G. J. Kistemaker, and B. J. Van Doormaal. 2000. Experience with a test-day model. J. Dairy Sci. 83:1135–1144.[Abstract]
Strandén, I., and D. Gianola. 1997. Gaussian versus Student’s-t mixed effects linear models for milk yield in Ayrshire Cattle. 48th Annual Meeting of the EAAP, Vienna, Austria.
Strandén, I., and D. Gianola. 1999. Mixed effects linear models with t-distributions for quantitative genetic analysis: A Bayesian approach. Genet. Sel. Evol. 31:25–34.
Wiggans, G. R., P. M. VanRaden, and J. C. Philpot. 2003. Detection and adjustment of abnormal test day yields. J. Dairy Sci. 86:2721–2724.[Abstract/Free Full Text]
TOP
ABSTRACT
INTRODUCTION
