Linear mixed-effects models and least confounded residuals in the modelling of Holstein calves’ performance

Some experimental studies are carried out considering a longitudinal feature. In view of this, the classical regression models are not able to handle with it, since the independence assumption between the observations is violated. To handle with this kind of data it was proposed the called linear mixed-effects models, where it is possible to model the response variable taking into account the correlation between the observations, and even between the response variables, when there are two or more of them in study, setting a bivariate or multivariate scenario, respectively. For the diagnosis of the linear mixed-effects models the least confounded residuals are quite recommended due to their lower bias in relation to other types of residuals, but it is not so used in the literature. Using a data set of dairy calves’ performance according to three different diets over eight weeks, the linear mixed-effects models theory under univariate and bivariate approach was applied alongside the least confounded residuals in the diagnosis of the model for both approaches. Comparing the univariate and bivariate approaches, the last one was more informative presenting lower standard errors’ values for its estimates, while the least confounded residuals was more efficient than the classical residuals present in the literature.


Introduction
The initial Holstein calves' performance is intrinsically related to several factors, such as the food composition (Berends et al., 2014).The liquid diet is the main nutrient source for calves; however, the solid diet presents a significant influence on the weight gain of the animals, which may lead to a productivity increase of those who consume great quantities of concentrate, mainly in the preweaning phase (Stamey et al., 2012).
Considering such scenario, the need of understanding the performance of dairy calves weight gain according to the solid diet is important to indicate a more adequate formulation to the animals.This need is more evident when taking into account the forager inclusion into the diet, since some studies indicate considerable benefits over the animal performance (Castells et al., 2012;Overvest et al., 2016;Bagheri et al., 2021).
An increase of the complexity at this kind of problem occurs when the interest lies on the study of the response variables over a determined time interval, considering the possible no null correlation between the observations taken repeatedly in the same individual.This context justifies the use of the univariate or multivariate linear mixed-effects models (Laird & Ware, 1982;Oskrochi et al., 2016).In the longitudinal experiments with animals, some disturbance in the data set can be expected, resultant of the loss of observations, causing an unbalancing on the data set, leading to a statistical efficiency loss and generating biased estimates (Zhu et al., 2014).
For the residuals diagnosis of the linear mixed-effects models, Hilden-Minton (1995) proposed the use of the least confounded residuals, which are not correlated to the confounded fraction associated to the parameters vector of random effects.Such residuals present lesser bias than the conditional and marginal ones, being highly recommended to the mixed models diagnosis (Loy & Hofmann, 2014).
The aim of this work was to study the weight and total solid feed intake behavior of the dairy calves submitted to three distinct diets over time, under univariate and bivariate approaches using linear mixed-effects models in a highly unbalanced data scenario.A secondary objective was to evaluate the use of the least confounded residuals in the fitted models diagnosis at both approaches when compared to the studentized conditional and marginal residuals.The tested hypotheses were: the bivariate approach better explains the longitudinal behavior of the two response variables and the use of the least confounded residuals is more adequate in the fitted models diagnosis.

Materials and Methods
The data set used in this work come from an experiment carried out in the Luiz de Queiroz College of Agriculture, University of São Paulo, Piracicaba -SP, that lasted 8 weeks.The objective of such experiment was studying the influence of the different non-forager fibre sources compared to the forager sources with low fibrous quality on the dairy calves' performance with several observed response variables (Poczynek et al., 2020).Thirty-five dairy calves distributed in twenty-two blocks, with the birth weight as the blocking factor, were used in the study.Three types of diets were compared: a diet containing 22% of neutral detergent fibre (NDF) without hay addition, a diet containing 22% of NDF with hay addition and diet containing 31% of NDF without hay addition.
Among the response variables evaluated in the experiment, the body weight (kg) and the weekly average total solid feed intake (g) of each animal were used at the present work .This because it is expected there is a non-null correlation between both.The animals were weekly weighted before the morning feeding and the solid feed intake measurements were taken daily before the new concentrate offer and calculated by the difference between the amount offered on the previous day and the leftover in the trough on time of the next offer.The three diets were randomized to the animals within each block.For being an animal inherent feature, it was studied the inclusion of the weight value in the random structure of the model overriding the block factor in the fixed and random structures of the model.
In the statistical analyses of the data the univariate and bivariate linear mixed-effects models were used.At the first stage and in both approaches, the maximal models were defined to the fixed and random structures.The polynomial degree that well explains the behavior of the response variable over time and the suggestion of an ideal random effects structure associated to its coefficients were embased in the visual analysis of the individual and average profiles.
In Figure 1 it is seen that the individual profiles of the calves weight present a curvy behavior that can be well explained by a second degree polynomial, as well as the average profile of each diet presented in Figure 3.Because of this aspect and the variability of the individual profiles, the maximum linear mixed model used to explain the longitudinal behavior of the calves weight submitted to three different diets is given by: where y ij is the response of the ith individual in the jth week considering the three diets (D  40, j = 1, 2, 3,..., 8; β 0 is the intercept of the curve related to the D 1 ; β 1 and β 2 are the first and second degree coefficients, respectively, of the polynomial related to the D 1 ; β 3 is the intercept of the curve related to the D 2 , while β 4 and β 5 are the first and second degree coefficients associated, respectively, to the same diet; β 6 is the intercept of the curve associated to the D 3 and β 7 and β 8 are the first and the second degree polynomial coefficients associated to the D 3 ; b 0 i is the random effect related to the intercept of the ith individual; b 1 i and b 2 i are the random effects related, respectively, to the first and second polynomial degree for the ith individual; ε ij is the random error associated to y ij .The individual and average profiles of solid feed intake presented in the Figure 2 and Figure 3, respectively, show a distinct behavior than that one presented for the calves weight.The solid feed intake observations through the weeks indicates the need of an up forth degree polynomial to explain such behavior.The proximity of the observations taken on the first week indicates that the inclusion of a distinct random effect to the intercept of the curves is not necessary.In view of this it was used the maximum model presented in (2), that considers a forth degree polynomial to explain the average behavior of the responses over the time of each diet and random effects related to all polynomial coefficients, except the intercept.
in this case β 0 is the coefficient of the intercept related to the D 1 , while β 1 , β 2 , β 3 and β 4 are up forth degree coefficients to the same diet.The meaning of the parameters associated to D 2 and D 3 is analog to the presented for D 1 .In the random structure, b 1 i and b 2 i have the same interpretation of those ones shown in (1) for the weight; b 3 1 and b 4 1 are, respectively, the random effects associated to the third and forth degree polynomial coefficients for the ith individual.
The maximum model proposed for the bivariate approach with both response variables is given by the expression (3): where y rij is the rth response of the ith individual in the jth week, being r 1 and r 2 the weight and the solid feed intake indexes, respectively and r only when the effect is shared by both responses; β 0r is the intercept related to the D 1 for the rth response, while β 1r and β 2r are the first and the second degree coefficients of the polynomial associated to the jth week for D 1 for the rth response; Associated to the D 2 effect, β 3 , β4 and β 5 are the intercept, first and second degree polynomial coefficients of the jth week for the rth response.The same meaning is given for the coefficients β 6 , β 7 and β 8 related to the D 3 ; β 9 and β 10 are the third and forth degree polynomial coefficients associated to the jth week only for the solid feed intake response.The same interpretation is given to β 11 and β 12 associated to the D 2 effect and β 13 and β 14 related to the D 3 effect; b 0r 1i is the random effect associated to the intercept for the ith individual weight; b 1r i and b 2r i are the random effects associated, respectively, to the first and the second degree polynomial coefficients for both responses on the ith individual; β 3r 2i and β 4r 2i are the random effects associated to the third and forth degree polynomial coefficients for the ith individual only for the solid feed intake; ε rij is the random error associated to y rij .In a second stage several random effects structures and different structures of inter-and intraindividuals covariance matrices were compared, aiming to choose the best set of random effects and the best covariance structure.On the comparison of the random effects structures it was used the called top-down strategy proposed by Verbeke et al. (1997).In all model comparisons the Likelihood-Ratio Test (Lehmann & Romano, 1986) and the Bayesian Information Criterion (BIC) (Schwarz, 1978) were used.
In a third stage, comparisons were made between the parameters of the fixed structure of the maximum model with simplest models according to evident suggestions by the coefficient estimates of the fitted polynomial for each diet and the fitted average curve plots.Admitting the possibility that a new random structure fits better than that one selected on the previous step, in a fourth stage, the select procedure for random effects was made using the same tests and statistics indicated on the second stage.
After the final choose of each univariate or bivariate model, its diagnoses were made using, mainly, residual plot analyses considering the least confounded residuals, studentized marginal and conditional residuals.
The marginal residuals are characterized as the difference between the vector of observed values and the estimated general mean, being used for the fixed effects linearity verification and the presence of outliers (Singer et al., 2013).According to the same authors, on the conditional residuals case, these consist as the difference between the vector of observed values and the predicted value by the model, getting the function of verifying the presence of outliers, as well as the normality and conditional errors homogeneity assumptions.However, both of them are not indicated for the linear models diagnoses, since they can present correlation between themselves, rupturing the independence assumption (Loy et al., 2013).In order to overcome this problem, such residuals should be standardized dividing them by its standard deviation estimates, and hereafter called as studentized residuals.Such procedure is done according shown in ( 4) and ( 5) expressions.
where r m and r c are the vectors of marginal and conditional residuals, respectively; r sc i and r sc i the studentized residuals; y the vector of observations; X and Z the design matrices for fixed and random factors, respectively; β the vector of estimated parameters associated to the fixed effects and b the vector of predicted random effects.Since it is two stages hierarchical model, the linear mixed model generates a vector of residuals that got a dependency with the vector of random effects, being necessary a more thorough diagnosis in view of such dependency (Loy et al., 2013).For that, Hilden-Minton (1995) proposes the called least confounded residuals that consists on an appropriate linear combination over the conditional residuals minimizing the confounding between the two stages of the model.Once such transformation is done, the diagnoses follow the usual protocols.
The Shapiro-Wilk (Shapiro & Wilk, 1965) and the Henze-Zirkler (Henze & Zirkler, 1990) tests were used to verify the normality adherence assumption of the residuals for univariate and bivariate approach, respectively.In both cases were also used the following plots: Q-Q plot to visualize the residuals normality; predicted values versus observed values dispersion, and plots with fitted average and individual curves.The study of the presence of atypical points in the bivariate analysis was embased in the bivariate boxplot (Rousseeuw et al., 1999).The comparison between the least confounded residuals, studentized marginal and conditional residuals was made embased in plots analyses in both approaches.The models were fitted in R language programming using the lme() function from the nlme package.5% of significance level was used for all hypothesis tests, which is quite common in studies carried out within agricultural sciences.

Results and Discussion
On the Table 1 is evident that the weight variance increases over time, as well as high and positive correlations between the measures taken in the different weeks.Furthermore, there is a decrease of the correlation as the interval between the weeks increases, which is expected in the longitudinal data.In an analog way to the weight, the Table 2 also enhance to the total solid feed intake an accentuated increase in the variability over time, lowest correlation coefficients and a decrease tendency in the sample correlation as the interval between the weeks increases.
Considering the bivariate approach, the Table 3 presents the correlation between the total solid feed intake and the weight of the animals over the weeks.As it is expected, in a general way the correlation tends to be greater when the two responses measures are taken in the same instant.It can also be realized that this feature is more evident in the last weeks than in the first ones, suggesting that over the weeks the calves developed a certain intake and weight standard.This also may be  The Table 4 presents the coefficient estimates and the parameters standard errors of the final fitted model to the calves weight under univariate approach that resulted, for the fixed structure, an only second degree curve, indicating that the average weight of the calves over the experimental period has the same behavior for both diets.The random structure included random effects for all polynomial coefficients.The best inter-individual covariance structure was the unstructured type, while the correlation and the intra-individual variability was modeled by the first order autorregressive with variance heterogeneity structure.The average and individual fitted curves presented in the Figure 4 evidence that the final linear mixed model managed to model the average as well as the individual weights over time very well.
The Figure 5 (a) presents the plots of the predicted values versus the residuals of the fitted model and the Figure 5 (b) the Q-Q plots for the three types of residuals.It realizes that the least confounded and the conditional residuals present a random distribution related to the predicted values.The Q-Q plots for both the last confounded residuals and the marginal residuals show an apparent normality adherence, what is confirmed by the Shapiro-Wilk test, with p-value equals to 0.360 and 0.713, respectively.However the conditional residuals presented a greater quantity of atypical points than the least confounded residuals and did not show adherence to the normal distribution as accused by the Shapiro-Wilk test, with p-value <0.0001.It should be noted that the studentized marginal residuals distribution regarding the predicted values evidences a strong tendency of linear increase and variability, confirming once again the bias which such residuals have.In view of the exposed, it can be concluded that the least confounded residuals shows a better performance in the diagnosis of the final model quality.
Still on the univariate approach, the Table 5 presents the coefficient estimates and standard errors of the parameters of the final model fitted to the total solid feed intake.For the fixed structure, the final model considered different intercepts for the diets containing 22% of NDF and another intercept for the diet with 31% of NDF.The same behavior occurred to the coefficient of second degree.The coefficients of first, third and fourth degrees were considered equals to the three diets.In the random structure, the random effect of the second degree of week was maintained.The chosen intra-individuals covariance structure was the autorregressive of first order with variance heterogeneity, such as the chosen one for the calve weights.The Figure 6 presents the fitted individual and average polynomial curves for total solid feed intake submitted to the two diets containing 22% of NDF and 31% of NDF.
Most of the individual curves fitted by the final model and presented in the Figure 6 (a), well describes the behavior of the observed responses for each calf.In the Figure 6 (b) it can be seen that the wavelike behavior of the average total solid feed intake was well explained by the final model.Similar results were found by Movahedi et al. (2017) in an experiment with Holstein calves that compared different diet compositions between concentrates and forages sources, indicating distinct behaviors of intake according to the diet provided to the animals.
In the Figure 7 (a) it can be seen that the least confounded residuals present a random distribution in relation to the predicted values; the marginal and conditional residuals present an apparent tendency or a great number of residuals superior to three, in absolute value, which indicate the presence of atypical points.In the Figure 7 (b) it is notable that the least confounded residuals meet the assumption of normality adherence, different of what occurs with the marginal and conditional residuals.The Shapiro-Wilk test for the least confounded residuals, marginal and conditional ones presented p-values equal to 0.507, <0.0001 and <0.0001, respectively, confirming such assumptions.
The bivariate study resulted in a model that considers for the average weight a single second degree polynomial curve to the three diets and for the average total solid feed intake a single fourth degree polynomial curve without the third degree coefficient, indicating that the animals that were fed with different diets presented the same behavior of weight and intake.The coefficient estimates of the final model are presented in the Table 6.
The random structure for the weight of the calves maintained the random effects associated to the intercept and the first degree coefficient.For the total solid feed intake, random effects associated to the second and forth degree coefficients were mantained.The chosen inter-individuals covariance structure was the unstructured and for the intra-individuals structure was the autorregressive of first order with variance heterogeneity, for both response variables.
It is understandable that the differences in the fixed and random structure of the final bivariate model in relation to the ones found in the univariate approaches to the weight and intake of the calves are justified by the correlation between the responses.
The ability of the bivariate model in using this dependency relationship produced a decrease in the values of standard errors of the polynomial coefficients that are common to the both approaches, as the case of the fitted polynomial to the animals weight.
The fitted individual and average curves to the responses in the eight weeks are present in the Figures 8 and 9, respectively.
In both figures it can be seen a well fit by the bivariate linear mixed model to the individual and average profiles of weight and solid feed intake.As the univariate approach, the final bivariate model better fitted to the weight data than the intake one.Regarding the response variables behavior, such results are according to the ones got by Engelking et al. (2020), who concluded that the inclusion or not inclusion of forages in the dairy calves diet in the pre-weaning period did not result significant differences in the weight and solid feed intake.7 is presented the correlation matrix between the observations of weight and solid feed intake over time, fitted by the final bivariate marginal model.It can be seen that such matrix shows a common behavior in longitudinal data resulting from experiments involving the animal performance, whereupon the non-null correlation values tendency to decrease as increase the interval between the measurements, such behavior is more evident after the fourth week of evaluation, corroborating with the descriptive results given in the Table 3.It also perceives median correlation values between the weight and the total solid feed intake in some weeks.These aspects reinforce the proper choice of random effects and inter-and intra-individuals covariance matrices that compound the chosen bivariate linear mixed model to explain the simultaneous longitudinal behavior for both responses.

Conclusions
The bivariate approach proved to be more appropriate than the univariate approach in the analysis of the weight and total solid feed intake data of calves.The simultaneous use of data of these two response variables caused alterations in the coefficient estimates of the polynomial model that explains the longitudinal behavior of the responses in its standard errors and also in the inferences made over the parameters.The least confounded residuals of the univariate and bivariate models presented lesser biased behavior than the studentized marginal and conditional residuals, as cited by the literature.

Figure 2 .
Figure 2. Individual profiles of the weekly solid feed intake of calves fed each evaluated diet.

Figure 3 .
Figure 3. Average profiles of the weekly body weight and solid feed intake of calves fed different diets.

Figure 5 .
Figure 5. (a) Plot of predicted values versus residuals and (b) QQ-plot for each kind of residuals generated by the final fitted for the weight of the calves submitted to different diets.

Figure 6 .
Figure 6.(a) Individual profiles and fitted individual of total solid feed intake and (b) Average profile and fitted average curve by the final model for the diets containing 22% of NDF and 31% of NDF.

Figure 7 .
Figure 7. (a) Plot of predicted values versus residuals and (b) QQ-plot for each kind of residuals generated by the final fitted model for the total solid feed intake on the calves submitted to different diets.

Table 6 .Figure 8 .Figure 9 .
Figure 8. Fitted individual and average of the weekly weight development under bivariate approach.

Figures 10 (
Figures 10 (a) and 10 (b) present the plots of predicted values versus the residuals generated by the bivariate model for the weight and solid feed intake, respectively, considering the least confounded, Individual profiles of the weekly weight of calves fed each evaluated diet.

Table 1 .
Variances (diagonal), covariances (under the diagonal) and correlation (above the diagonal) of the calves weight over eight weeks

Table 2 .
Variances (diagonal), covariances (under the diagonal) and correlation (above the diagonal) of the calves total solid feed average intake between the weeks explained by the low solid feed intake during the first weeks of age, when calves rely mostly on the liquid diet for growth.

Table 3 .
Correlation between the calves total solid feed average intake and the weight over the weeks

Table 5 .
Coefficient estimates and parameter standard errors of the second degree polynomial that explains the calves weight behavior

Table 7 .
Bivariate correlation matrix for the marginal model