Statistical analysis of test-day milk yields using random regression models for the comparison of feeding groups during the lactation period

Abstract

Random regression models are widely used in the field of animal breeding for the genetic evaluation of daily milk yields from different test days. These models are capable of handling different environmental effects on the respective test day, and they describe the characteristics of the course of the lactation period by using suitable covariates with fixed and random regression coefficients. As the numerically expensive estimation of parameters is already part of advanced computer software, modifications of random regression models will considerably grow in importance for statistical evaluations of nutrition and behaviour experiments with animals. Random regression models belong to the large class of linear mixed models. Thus, when choosing a model, or more precisely, when selecting a suitable covariance structure of the random effects, the information criteria of Akaike and Schwarz can be used. In this study, the fitting of random regression models for a statistical analysis of a feeding experiment with dairy cows is illustrated under application of the program package SAS. For each of the feeding groups, lactation curves modelled by covariates with fixed regression coefficients are estimated simultaneously. With the help of the fixed regression coefficients, differences between the groups are estimated and then tested for significance. The covariance structure of the random and subject-specific effects and the serial correlation matrix are selected by using information criteria and by estimating correlations between repeated measurements. For the verification of the selected model and the alternative models, mean values and standard deviations estimated with ordinary least square residuals are used.     

Application of the covariance function approach with an iterative two-stage algorithm to the estimation of parameters of a random regression test day model for dairy production traits

The covariance function approach with an iterative two-stage algorithm of LIU et al. (2000) was applied to estimate parameters for the Polish Black-and-White dairy population based on a sample of 338 808 test day records for milk, fat, and protein yields. A multiple trait sire model was used to estimate covariances of lactation stages. A third-order Legendre polynomial was subsequently fitted to the estimated (co)variances to derive (co)variances of random regression coefficients for both additive genetic and permanent environment effects. Daily and 305-day heritability estimates obtained are consistent with several studies which used both fixed and random regression test day models. Genetic correlations between any two days in milk (DIM) of the same lactation as well as genetic correlations between the same DIM of two lactations were within a biologically acceptable range. It was shown that the applied estimation procedure can utilise very large data sets and give plausible estimates of (co)variance components.     

Analysis of water intake, dry matter intake and daily milk yield using different error covariance structures

The aim of the present study was to investigate the daily measured traits milk yield, water intake and dry matter intake with fixed and random regression models added with different error covariance structures. It was analysed whether these models deliver better model fitting in contrast to conventional fixed and random regression models. Furthermore, possible autocorrelation between repeated measures was investigated. The effect of model choice on statistical inference was also tested. Data recording was performed on the Futterkamp dairy research farm of the Chamber of Agriculture of Schleswig-Holstein. A dataset of about 21 000 observations from 178 Holstein cows was used. Average milk yield, water intake and dry matter intake were 34.9, 82.4 and 19.8 kg, respectively. Statistical analysis was performed using different linear mixed models. Lactation number, test day and the parameters to model the function of lactation day were included as fixed effects. Different structures were tested for the residuals; they were compared for their ability to fit the model using the likelihood ratio test, and Akaike's and Bayesian's information criteria. Different autocorrelation patterns were found. Adjacent repeated measures of daily milk yield were highest correlated (p1 = 0.32) in contrast to measures further apart, while for water intake and dry matter intake, the measurements with a lag of two units had the highest correlations with p2 = 0.11 and 0.12. The covariance structure of TOEPLITZ was most suitable to indicate the dependencies of the repeated measures for all traits. Generally, the most complex model, random regression with the additional covariance structure TOEPLITZ(4), provided the lowest information criteria. Furthermore, the model choice influenced the significance values of one fixed effect and therefore the general inference of the data analysis. Thus, the random regression + TOEPLITZ(4) model is recommended for use for the analysis of equally spaced datasets of milk yield, water intake and dry matter intake.     

Analysis of multi-type recurrent events in longitudinal studies; application to a skin cancer prevention trial

We consider the statistical modeling and analysis of replicated multi-type point process data with covariates. Such data arise when heterogeneous subjects experience repeated events or failures which may be of several distinct types. The underlying processes are modeled as nonhomogeneous mixed Poisson processes with random (subject) and fixed (covariate) effects. The method of maximum likelihood is used to obtain estimates and standard errors of the failure rate parameters and regression coefficients. Score tests and likelihood ratio statistics are used for covariate selection. A graphical test of goodness of fit of the selected model is based on generalized residuals. Measures for determining the influence of an individual observation on the estimated regression coefficients and on the score test statistic are developed. An application is described to a large ongoing randomized controlled clinical trial for the efficacy of nutritional supplements of selenium for the prevention of two types of skin cancer.     

Random regression models using different functions to model test-day milk yield of Brazilian Holstein cows

We analyzed 152,145 test-day records from 7317 first lactations of Holstein cows recorded from 1995 to 2003. Our objective was to model variations in test-day milk yield during the first lactation of Holstein cows by random regression model (RRM), using various functions in order to obtain adequate and parsimonious models for the estimation of genetic parameters. Test-day milk yields were grouped into weekly classes of days in milk, ranging from 1 to 44 weeks. The contemporary groups were defined as herd-test-day. The analyses were performed using a single-trait RRM, including the direct additive, permanent environmental and residual random effects. In addition, contemporary group and linear and quadratic effects of the age of cow at calving were included as fixed effects. The mean trend of milk yield was modeled with a fourth-order orthogonal Legendre polynomial. The additive genetic and permanent environmental covariance functions were estimated by random regression on two parametric functions, Ali and Schaeffer and Wilmink, and on B-spline functions of days in milk. The covariance components and the genetic parameters were estimated by the restricted maximum likelihood method. Results from RRM parametric and B-spline functions were compared to RRM on Legendre polynomials and with a multi-trait analysis, using the same data set. Heritability estimates presented similar trends during mid-lactation (13 to 31 weeks) and between week 37 and the end of lactation, for all RRM. Heritabilities obtained by multi-trait analysis were of a lower magnitude than those estimated by RRM. The RRMs with a higher number of parameters were more useful to describe the genetic variation of test-day milk yield throughout the lactation. RRM using B-spline and Legendre polynomials as base functions appears to be the most adequate to describe the covariance structure of the data.     

Correcting for measurement error in individual-level covariates in nonlinear mixed effects models

The nonlinear mixed effects model is used to represent data in pharmacokinetics, viral dynamics, and other areas where an objective is to elucidate associations among individual-specific model parameters and covariates; however, covariates may be measured with error. For additive measurement error, we show substitution of mismeasured covariates for true covariates may lead to biased estimators for fixed effects and random effects covariance parameters, while regression calibration may eliminate bias in fixed effects but fail to correct that in covariance parameters. We develop methods to take account of measurement error that correct this bias and may be implemented with standard software, and we demonstrate their utility via simulation and application to data from a study of HIV dynamics.     

Effects of residual smoothing on the posterior of the fixed effects in disease-mapping models

Disease-mapping models for areal data often have fixed effects to measure the effect of spatially varying covariates and random effects with a conditionally autoregressive (CAR) prior to account for spatial clustering. In such spatial regressions, the objective may be to estimate the fixed effects while accounting for the spatial correlation. But adding the CAR random effects can cause large changes in the posterior mean and variance of fixed effects compared to the nonspatial regression model. This article explores the impact of adding spatial random effects on fixed effect estimates and posterior variance. Diagnostics are proposed to measure posterior variance inflation from collinearity between the fixed effect covariates and the CAR random effects and to measure each region's influence on the change in the fixed effect's estimates by adding the CAR random effects. A new model that alleviates the collinearity between the fixed effect covariates and the CAR random effects is developed and extensions of these methods to point-referenced data models are discussed.     

多反应变量重复测量的混合效应模型及其应用

本文给出了多反应变量重复测量的协方差矩阵结构,探讨了用迭代广义最小二乘法来求解其带协变量和不带协变量的混合效应模型中固定效应和随机效应系数,并对１９９１年四川省高血压调查资料进行实例分析,得到其结论符合实际情况．     

Data transformation for rank reduction in multi-trait MACE model for international bull comparison

Since many countries use multiple lactation random regression test day models in national evaluations for milk production traits, a random regression multiple across-country evaluation (MACE) model permitting a variable number of correlated traits per country should be used in international dairy evaluations. In order to reduce the number of within country traits for international comparison, three different MACE models were implemented based on German daughter yield deviation data and compared to the random regression MACE. The multiple lactation MACE model analysed daughter yield deviations on a lactation basis reducing the rank from nine random regression coefficients to three lactations. The lactation breeding values were very accurate for old bulls, but not for the youngest bulls with daughters with short lactations. The other two models applied principal component analysis as the dimension reduction technique: one based on eigenvalues of a genetic correlation matrix and the other on eigenvalues of a combined lactation matrix. The first one showed that German data can be transformed from nine traits to five eigenfunctions without losing much accuracy in any of the estimated random regression coefficients. The second one allowed performing rank reductions to three eigenfunctions without having the problem of young bulls with daughters with short lactations.     

Isolating the cow-specific part of residual energy intake in lactating dairy cows using random regressions

The ability to properly assess and accurately phenotype true differences in feed efficiency among dairy cows is key to the development of breeding programs for improving feed efficiency. The variability among individuals in feed efficiency is commonly characterised by the residual intake approach. Residual feed intake is represented by the residuals of a linear regression of intake on the corresponding quantities of the biological functions that consume (or release) energy. However, the residuals include both, model fitting and measurement errors as well as any variability in cow efficiency. The objective of this study was to isolate the individual animal variability in feed efficiency from the residual component. Two separate models were fitted, in one the standard residual energy intake (REI) was calculated as the residual of a multiple linear regression of lactation average net energy intake (NEI) on lactation average milk energy output, average metabolic BW, as well as lactation loss and gain of body condition score. In the other, a linear mixed model was used to simultaneously fit fixed linear regressions and random cow levels on the biological traits and intercept using fortnight repeated measures for the variables. This method split the predicted NEI in two parts: one quantifying the population mean intercept and coefficients, and one quantifying cow-specific deviations in the intercept and coefficients. The cow-specific part of predicted NEI was assumed to isolate true differences in feed efficiency among cows. NEI and associated energy expenditure phenotypes were available for the first 17 fortnights of lactation from 119 Holstein cows; all fed a constant energy-rich diet. Mixed models fitting cow-specific intercept and coefficients to different combinations of the aforementioned energy expenditure traits, calculated on a fortnightly basis, were compared. The variance of REI estimated with the lactation average model represented only 8% of the variance of measured NEI. Among all compared mixed models, the variance of the cow-specific part of predicted NEI represented between 53% and 59% of the variance of REI estimated from the lactation average model or between 4% and 5% of the variance of measured NEI. The remaining 41% to 47% of the variance of REI estimated with the lactation average model may therefore reflect model fitting errors or measurement errors. In conclusion, the use of a mixed model framework with cow-specific random regressions seems to be a promising method to isolate the cow-specific component of REI in dairy cows.     

Variable selection – A review and recommendations for the practicing statistician

Statistical models support medical research by facilitating individualized outcome prognostication conditional on independent variables or by estimating effects of risk factors adjusted for covariates. Theory of statistical models is well‐established if the set of independent variables to consider is fixed and small. Hence, we can assume that effect estimates are unbiased and the usual methods for confidence interval estimation are valid. In routine work, however, it is not known a priori which covariates should be included in a model, and often we are confronted with the number of candidate variables in the range 10–30. This number is often too large to be considered in a statistical model. We provide an overview of various available variable selection methods that are based on significance or information criteria, penalized likelihood, the change‐in‐estimate criterion, background knowledge, or combinations thereof. These methods were usually developed in the context of a linear regression model and then transferred to more generalized linear models or models for censored survival data. Variable selection, in particular if used in explanatory modeling where effect estimates are of central interest, can compromise stability of a final model, unbiasedness of regression coefficients, and validity of p‐values or confidence intervals. Therefore, we give pragmatic recommendations for the practicing statistician on application of variable selection methods in general (low‐dimensional) modeling problems and on performing stability investigations and inference. We also propose some quantities based on resampling the entire variable selection process to be routinely reported by software packages offering automated variable selection algorithms.     

The effect of retrospective sampling on binary regression models for clustered data

Recently a great deal of attention has been given to binary regression models for clustered or correlated observations. The data of interest are of the form of a binary dependent or response variable, together with independent variables X1,...., Xk, where sets of observations are grouped together into clusters. A number of models and methods of analysis have been suggested to study such data. Many of these are extensions in some way of the familiar logistic regression model for binary data that are not grouped (i.e., each cluster is of size 1). In general, the analyses of these clustered data models proceed by assuming that the observed clusters are a simple random sample of clusters selected from a population of clusters. In this paper, we consider the application of these procedures to the case where the clusters are selected randomly in a manner that depends on the pattern of responses in the cluster. For example, we show that ignoring the retrospective nature of the sample design, by fitting standard logistic regression models for clustered binary data, may result in misleading estimates of the effects of covariates and the precision of estimated regression coefficients.     

Continuous time Markov models for binary longitudinal data

Longitudinal data usually consist of a number of short time series. A group of subjects or groups of subjects are followed over time and observations are often taken at unequally spaced time points, and may be at different times for different subjects. When the errors and random effects are Gaussian, the likelihood of these unbalanced linear mixed models can be directly calculated, and nonlinear optimization used to obtain maximum likelihood estimates of the fixed regression coefficients and parameters in the variance components. For binary longitudinal data, a two state, non-homogeneous continuous time Markov process approach is used to model serial correlation within subjects. Formulating the model as a continuous time Markov process allows the observations to be equally or unequally spaced. Fixed and time varying covariates can be included in the model, and the continuous time model allows the estimation of the odds ratio for an exposure variable based on the steady state distribution. Exact likelihoods can be calculated. The initial probability distribution on the first observation on each subject is estimated using logistic regression that can involve covariates, and this estimation is embedded in the overall estimation. These models are applied to an intervention study designed to reduce children's sun exposure.     

Breeding value accuracy estimates for growth traits using random regression and multi-trait models in Nelore cattle

We quantified the potential increase in accuracy of expected breeding value for weights of Nelore cattle, from birth to mature age, using multi-trait and random regression models on Legendre polynomials and B-spline functions. A total of 87,712 weight records from 8144 females were used, recorded every three months from birth to mature age from the Nelore Brazil Program. For random regression analyses, all female weight records from birth to eight years of age (data set I) were considered. From this general data set, a subset was created (data set II), which included only nine weight records: at birth, weaning, 365 and 550 days of age, and 2, 3, 4, 5, and 6 years of age. Data set II was analyzed using random regression and multi-trait models. The model of analysis included the contemporary group as fixed effects and age of dam as a linear and quadratic covariable. In the random regression analyses, average growth trends were modeled using a cubic regression on orthogonal polynomials of age. Residual variances were modeled by a step function with five classes. Legendre polynomials of fourth and sixth order were utilized to model the direct genetic and animal permanent environmental effects, respectively, while third-order Legendre polynomials were considered for maternal genetic and maternal permanent environmental effects. Quadratic polynomials were applied to model all random effects in random regression models on B-spline functions. Direct genetic and animal permanent environmental effects were modeled using three segments or five coefficients, and genetic maternal and maternal permanent environmental effects were modeled with one segment or three coefficients in the random regression models on B-spline functions. For both data sets (I and II), animals ranked differently according to expected breeding value obtained by random regression or multi-trait models. With random regression models, the highest gains in accuracy were obtained at ages with a low number of weight records. The results indicate that random regression models provide more accurate expected breeding values than the traditionally finite multi-trait models. Thus, higher genetic responses are expected for beef cattle growth traits by replacing a multi-trait model with random regression models for genetic evaluation. B-spline functions could be applied as an alternative to Legendre polynomials to model covariance functions for weights from birth to mature age.     

Quantifying Spatial Disparities in Neonatal Mortality Using a Structured Additive Regression Model

Background

Neonatal mortality contributes a large proportion towards early childhood mortality in developing countries, with considerable geographical variation at small areas within countries.

Methods

A geo-additive logistic regression model is proposed for quantifying small-scale geographical variation in neonatal mortality, and to estimate risk factors of neonatal mortality. Random effects are introduced to capture spatial correlation and heterogeneity. The spatial correlation can be modelled using the Markov random fields (MRF) when data is aggregated, while the two dimensional P-splines apply when exact locations are available, whereas the unstructured spatial effects are assigned an independent Gaussian prior. Socio-economic and bio-demographic factors which may affect the risk of neonatal mortality are simultaneously estimated as fixed effects and as nonlinear effects for continuous covariates. The smooth effects of continuous covariates are modelled by second-order random walk priors. Modelling and inference use the empirical Bayesian approach via penalized likelihood technique. The methodology is applied to analyse the likelihood of neonatal deaths, using data from the 2000 Malawi demographic and health survey. The spatial effects are quantified through MRF and two dimensional P-splines priors.

Results

Findings indicate that both fixed and spatial effects are associated with neonatal mortality.

Conclusions

Our study, therefore, suggests that the challenge to reduce neonatal mortality goes beyond addressing individual factors, but also require to understanding unmeasured covariates for potential effective interventions.     

Identifying multispecies synchrony in response to environmental covariates

The importance of multispecies models for understanding complex ecological processes and interactions is beginning to be realized. Recent developments, such as those by Lahoz‐Monfort et al. (2011), have enabled synchrony in demographic parameters across multiple species to be explored. Species in a similar environment would be expected to be subject to similar exogenous factors, although their response to each of these factors may be quite different. The ability to group species together according to how they respond to a particular measured covariate may be of particular interest to ecologists. We fit a multispecies model to two sets of similar species of garden bird monitored under the British Trust for Ornithology's Garden Bird Feeding Survey. Posterior model probabilities were estimated using the reversible jump algorithm to compare posterior support for competing models with different species sharing different subsets of regression coefficients. There was frequently good agreement between species with small asynchronous random‐effect components and those with posterior support for models with shared regression coefficients; however, this was not always the case. When groups of species were less correlated, greater uncertainty was found in whether regression coefficients should be shared or not. The methods outlined in this study can test additional hypotheses about the similarities or synchrony across multiple species that share the same environment. Through the use of posterior model probabilities, estimated using the reversible jump algorithm, we can detect multispecies responses in relation to measured covariates across any combination of species and covariates under consideration. The method can account for synchrony across species in relation to measured covariates, as well as unexplained variation accounted for using random effects. For more flexible, multiparameter distributions, the support for species‐specific parameters can also be measured.     

Genetic relationship between milk urea nitrogen and reproductive performance in Holstein dairy cows

The objective of this study was to describe the genetic and phenotypic relationship between milk urea nitrogen (MUN) and reproductive traits in Iranian Holstein dairy cows. Test-day MUN data obtained from 57 301 dairy cows on 20 large dairy herds in Iran between January 2005 and June 2009. Genetic parameters for MUN and reproductive traits were estimated with a five-trait model using ASREML program. Random regression test-day models were used to estimate heritabilities separately for MUN from first, second and third lactations. Regression curves were modeled using Legendre polynomials of order 3. Herd-year-season along with age at calving was included as fixed effects in all models for reproductive traits. Heritabilities for MUN and reproductive traits were estimated separately for first lactation, second lactation and third lactation. The estimated heritabilities for MUN varied from 0.18 to 0.22. The heritability estimate was low for reproductive traits, which ranged from 0.02 to 0.06 for different traits and across parities. Except for days open, phenotypic and genetic correlations of MUN with reproductive performance traits were close to zero. Genetic correlations between MUN and days open were 0.23, 0.35 and 0.45 in first, second and third lactation, respectively. However, the phenotypic correlation between MUN at different parities was moderate (0.28 to 0.35), but the genetic correlation between MUN at different parities was high and ranged from 0.84 to 0.97. This study shows a limited application of MUN for use in selection programs to improve reproductive performance.     

Multiple trait model combining random regressions for daily feed intake with single measured performance traits of growing pigs

A random regression model for daily feed intake and a conventional multiple trait animal model for the four traits average daily gain on test (ADG), feed conversion ratio (FCR), carcass lean content and meat quality index were combined to analyse data from 1 449 castrated male Large White pigs performance tested in two French central testing stations in 1997. Group housed pigs fed ad libitum with electronic feed dispensers were tested from 35 to 100 kg live body weight. A quadratic polynomial in days on test was used as a regression function for weekly means of daily feed intake and to escribe its residual variance. The same fixed (batch) and random (additive genetic, pen and individual permanent environmental) effects were used for regression coefficients of feed intake and single measured traits. Variance components were estimated by means of a Bayesian analysis using Gibbs sampling. Four Gibbs chains were run for 550 000 rounds each, from which 50 000 rounds were discarded from the burn-in period. Estimates of posterior means of covariance matrices were calculated from the remaining two million samples. Low heritabilities of linear and quadratic regression coefficients and their unfavourable genetic correlations with other performance traits reveal that altering the shape of the feed intake curve by direct or indirect selection is difficult.     

Assessing variation in life‐history tactics within a population using mixture regression models: a practical guide for evolutionary ecologists

Mixed models are now well‐established methods in ecology and evolution because they allow accounting for and quantifying within‐ and between‐individual variation. However, the required normal distribution of the random effects can often be violated by the presence of clusters among subjects, which leads to multi‐modal distributions. In such cases, using what is known as mixture regression models might offer a more appropriate approach. These models are widely used in psychology, sociology, and medicine to describe the diversity of trajectories occurring within a population over time (e.g. psychological development, growth). In ecology and evolution, however, these models are seldom used even though understanding changes in individual trajectories is an active area of research in life‐history studies. Our aim is to demonstrate the value of using mixture models to describe variation in individual life‐history tactics within a population, and hence to promote the use of these models by ecologists and evolutionary ecologists. We first ran a set of simulations to determine whether and when a mixture model allows teasing apart latent clustering, and to contrast the precision and accuracy of estimates obtained from mixture models versus mixed models under a wide range of ecological contexts. We then used empirical data from long‐term studies of large mammals to illustrate the potential of using mixture models for assessing within‐population variation in life‐history tactics. Mixture models performed well in most cases, except for variables following a Bernoulli distribution and when sample size was small. The four selection criteria we evaluated [Akaike information criterion (AIC), Bayesian information criterion (BIC), and two bootstrap methods] performed similarly well, selecting the right number of clusters in most ecological situations. We then showed that the normality of random effects implicitly assumed by evolutionary ecologists when using mixed models was often violated in life‐history data. Mixed models were quite robust to this violation in the sense that fixed effects were unbiased at the population level. However, fixed effects at the cluster level and random effects were better estimated using mixture models. Our empirical analyses demonstrated that using mixture models facilitates the identification of the diversity of growth and reproductive tactics occurring within a population. Therefore, using this modelling framework allows testing for the presence of clusters and, when clusters occur, provides reliable estimates of fixed and random effects for each cluster of the population. In the presence or expectation of clusters, using mixture models offers a suitable extension of mixed models, particularly when evolutionary ecologists aim at identifying how ecological and evolutionary processes change within a population. Mixture regression models therefore provide a valuable addition to the statistical toolbox of evolutionary ecologists. As these models are complex and have their own limitations, we provide recommendations to guide future users.     

Comparison of One‐Step and Two‐Step Meta‐Analysis Models Using Individual Patient Data

The problem of combining information from separate trials is a key consideration when performing a meta‐analysis or planning a multicentre trial. Although there is a considerable journal literature on meta‐analysis based on individual patient data (IPD), i.e. a one‐step IPD meta‐analysis, versus analysis based on summary data, i.e. a two‐step IPD meta‐analysis, recent articles in the medical literature indicate that there is still confusion and uncertainty as to the validity of an analysis based on aggregate data. In this study, we address one of the central statistical issues by considering the estimation of a linear function of the mean, based on linear models for summary data and for IPD. The summary data from a trial is assumed to comprise the best linear unbiased estimator, or maximum likelihood estimator of the parameter, along with its covariance matrix. The setup, which allows for the presence of random effects and covariates in the model, is quite general and includes many of the commonly employed models, for example, linear models with fixed treatment effects and fixed or random trial effects. For this general model, we derive a condition under which the one‐step and two‐step IPD meta‐analysis estimators coincide, extending earlier work considerably. The implications of this result for the specific models mentioned above are illustrated in detail, both theoretically and in terms of two real data sets, and the roles of balance and heterogeneity are highlighted. Our analysis also shows that when covariates are present, which is typically the case, the two estimators coincide only under extra simplifying assumptions, which are somewhat unrealistic in practice.