Animals and data set descriptions

The data used in the current study were collected from the goat herd of the Goat Sector at the Federal University of Viçosa (Universidade Federal de Viçosa-UFV), Viçosa, MG, Brazil. A total of 14 003 weighing records from female goats were collected between 1992 and 2010. Only those goats with four or more recorded weight measurements were included. Data from animals without recorded birth dates were excluded. Then the data were divided according to the genotype of the animals: Saanen (+S) and Alpine (+A). A graphical analysis was conducted to assess the consistency of the data and to identify the excluded end points, including data from animals of all other genotypes. The database is summarized in Table 1.

Table 1.

Summary of the data set assessed

Animal management consisted of a free-stall milk production system in which animals had access to maize silage and concentrate ad libitum. There are normally two calving seasons per year for goats; however, each animal was allowed just one calving per year. During the growth phase, the animals were weaned, on average, when they were 18 kg (110–140 days old) and the first pregnancy occurred, on average, at 35 kg body weight (340–370 days old). The goats lactated for c. 300 days, and their average milk production was 2·6 kg milk/day. The number of lactations ranged from zero to seven.

The number of weighings by animal ranged from 4 to 21 (Table 1) and the minimum and maximum differences between two consecutive weighings were 1 and 1190 days, respectively. These features make the data set extremely unbalanced, spatially and numerically. Furthermore, the initial graph analysis suggested that there was an increase in the variance of the data over time (the independent parameter). Thus, a special approach to modelling the growth curves was required.

Growth curve functions

Five nonlinear models commonly used in the literature for describing animal growth curves were evaluated. In these models, shown in Table 2, y is the animal weight at time t _i; β ₁ the estimated body weight of the animal when t→∞ or asymptotic body weight; β ₂ the constant of integration; β ₃ the rate constant, which determines the spread of the curve along the time axis; and β ₄ is the constant that determines the proportion of the final size at which the inflection point (IP) occurs. The IP can be determined using the equation $IP_{(\% \beta _1 )} = \beta _4 ^{1/(1 - \beta _4 )} $, and the instantaneous absolute growth rate can be calculated using the equation $IAGR = \beta _3 \times Y_{{\rm est}} [(\beta _1 /Y_{{\rm est}} )^{1 - \beta _4} - 1]/(1 - \beta _4 )$, where Y _est is the body weight estimated at a given time (Richards Reference Richards1959).

Table 2.

Evaluated growth curve functions

Nonlinear mixed model

A nonlinear mixed model methodology was used to obtain the growth curve parameters according to the following expression:

(1)$${\bf y}_i = f\,({\bf x}_{ij} {,}{\bi \beta}, {\bi u}_i ) + {\bi e}_i $$

where y_i is a vector (n _i×1) of recorded body weights for the subject i=1, 2,…, m; m is the number of goats; f is a nonlinear function of the covariate matrix x_ij; x_ij is a matrix (n _i×2) of independent variables, in which the first column contains the ith age of measurement and the second column contains the jth subject (goat) whose body weight was measured; β is a vector (p _i×1) of unknown fixed-effect parameters; u _i is a vector (q _i×1) of unknown random-effect parameters unique to the subject goat i and assumed to follow a multivariate distribution with a mean of zero and a conditional unstructured covariance matrix, n×n (G); and e_i is a vector (n _i×1) of the error term assumed to follow a multivariate random normal distribution with zero mean and conditional covariance matrix structure (R).

The dependent variable (y) in this case is influenced by quantitative (time) and qualitative (genotype) variables. The effects of the quantitative variables are analysed via nonlinear regression. However, it is also necessary to assess the effects of the qualitative variable on the dependent variable and the regression parameters; thus, dummy variables z ₁ and z ₂ were created, where for the +S genotype, z ₁=0 and z ₂=1, and for the +A genotype, z ₁=1 and z ₂=0. With this approach, it was possible to estimate growth curve parameters (β ₁, β ₂, β ₃ and β ₄) independently for each genotype.

The growth curve modelling consisted of two steps, as described below. First, a set of candidate nonlinear mixed models were formulated that were preliminarily analysed to choose the function with the best fit. This first step was used to evaluate which curve parameters should have a random effect component and compared them with the traditional fixed-effect model. However, despite the high hardiness with which nonlinear mixed model methodology treats the correlated errors, it is important to note that heterogeneity and correlated errors can occur even with the inclusion of random effects (Meng & Huang Reference Meng and Huang2010; Yang & Huang Reference Yang and Huang2011). Thus, a refinement of the model chosen in the first step was necessary to accommodate the violated model assumptions (second step).

First step

A central phase in the model-building of mixed-effects models is to decide which of the parameters in the nonlinear model require random effects to account for their between-subject variations and which can be treated as purely fixed effects (Pinheiro & Bates Reference Pinheiro and Bates2000). Thus, the random parameters u ₁, u ₂ and u ₃ were created to account for between-subject variations for the β ₁, β ₂ and β ₃ parameters, respectively, enabling curve parameters to be estimated for each goat (Littell et al. Reference Littell, Milliken, Stroup, Wolfinger and Schabenberger2006). No random parameter was created for the β ₄ parameter in the Richards function because it was difficult to achieve the convergence criterion and the high correlations between the estimated parameters (overparameterized model). The following four combinations of random effects were evaluated: all of the parameters related to random effects (u ₁, u ₂ and u ₃), only parameters related to β ₁ and β ₂ (u ₁ and u ₂), only parameters related to β ₁ and β ₃ (u ₁ and u ₃), and only that related to β ₁ (u ₁). To avoid problems with floating-point errors and overflows, β ₃ was rescaled (×1000) in the models due to the difference in magnitude of this parameter in relation to the others.

First, the equations were fitted to each model and previously described combinations of random effects through PROC NLMIXED (SAS 9.3, see Appendix). The FIRO method was used to achieve convergence. The FIRO method uses the First-Order Method described by Beal & Sheiner (Reference Beal and Sheiner1982) to compute the integral over the random effects to the marginal maximum likelihood. The fixed-effects model was also estimated using PROC NLMIXED, except that, in this case, variance components were not included. The FIRO method does not support the model without random effects; thus, the Adaptive Gaussian Quadrature described by Pinheiro & Bates (Reference Pinheiro and Bates1995) was used. The estimated parameters, standard error and residual variance for the fixed-effects model were similar to those of the model obtained using PROC NLIN.

The following criteria were adopted to select the function that best described the growth curve of the goats: (1) convergence (the interactive process in NLMIXED converges at the nth interaction when $(g_k H_k ^{ - 1} g_k )/(|l_k | + 1\,e^{ - 6} ) \lt 1\,e^{ - 8} )$(SAS Institute Inc. 2008)); (2) the final Hessian matrix is non-singular and positive definite; (3) the approach described by Burnham & Anderson (Reference Burnham and Anderson2002) and Vieira et al. (Reference Vieira, Campos, da Silva, Tedeschi and Tamy2012) using the Akaike information criterion (AIC), the difference among AIC values (Δ_r), the Akaike weights or likelihood probabilities (w _r) and the evidence ratio or relative likelihood (ER_r), which can be computed using the following equations:

$${\rm AIC} = 2f\,({ \widehat\theta} ) + 2p$$

$$\Delta _{\rm r} = {\rm AIC}_{\rm r} - \min {\rm AIC}_{\rm r} $$

$$w_{\rm r} = \displaystyle{{\exp ( - \Delta _{\rm r} /2)} \over {\sum\nolimits_{r = 1}^R {\exp ( - \Delta _{\rm r} /2)}}} $$

$${\rm ER}_{\rm r} = \displaystyle{{{\rm max} w_{\rm r}} \over {w_{\rm r}}} $$

where f() is the negative of the marginal log-likelihood function, $\overset{ {\widehat \theta} $ is the vector of parameter estimate, and p is the number of parameters (SAS Institute Inc. 2008); and (4) a graphical analysis of the Pearson residuals against the predicted values was used to evaluate the model assumptions. The Pearson residuals were obtained for the fixed model and random model without a correlation structure as $r_i = e_i /\sqrt {\widehat{\sigma _e^2}} $ and for the model with a correlation structure (second step described below) as $r_i = \hat C^{' - 1} (y_i - \hat y_i )$, where $\hat C$ denotes the Cholesky root of the estimated R matrix (SAS Institute Inc. 2008).

Second step

The model with the best goodness-of-fit (combination between nonlinear equation and random parameters) was selected for further analysis.

Despite there being no sign of heteroscedasticity in the residuals of the chosen model, it was easy to identify high Pearson residuals obtained even when using random effects in the model. The higher values found for the Pearson residuals, many of them above three, were an indication of problems in the model. The first remedy would be to exclude the data from the data set and reparameterize the model, which would solve the problem but cause a loss of information.

Therefore, it was decided to use two tools to address this problem and the remaining correlated errors. The residual variance was modelled using the power-of-the-mean variance (PV) function. Using the notation described by Littell et al. (Reference Littell, Milliken, Stroup, Wolfinger and Schabenberger2006), the variance matrix R is assumed to be of the form ${\bi R} = diag(\sigma ^2 |x_i ^' \beta ^ * |^\theta )$, where θ is the power to be estimated. In order to choose the value of θ, the power that minimized the approximate −2log-likelihood was calculated (Littell et al. Reference Littell, Milliken, Stroup, Wolfinger and Schabenberger2006); this θ value was estimated to be 2·91.

As previously noted, correlated errors can occur even when random effects are used in the model. These errors primarily occur if there are parameters in the function that do not have a random parameter to account for their between-subject variability. An error structure matrix can be modelled to account for this interrelation. When the model is fitted without modelling the serial correlation, the independently and identically distributed error structure is assumed to be Var[e _i]=σ _e² and cov[e _i,e_j]=0. However, there is one intrinsic characteristic of the current data set that must be considered, namely, the unequally spaced measurement intervals for the animal subjects.

There is a range of structures available in the literature to address this issue. Structures that accept unequally spaced data have been described by Littell et al. (Reference Littell, Milliken, Stroup, Wolfinger and Schabenberger2006) and the SAS Institute Inc. (2008). Normally, the covariance is assumed to be a function of the distance between locations. If d _ij denotes the distance between the measurements made on the same animal, the covariance models have the general form COV[e _i,e_j]=σ_e²[f(d _ij)] (Littell et al. Reference Littell, Milliken, Stroup, Wolfinger and Schabenberger2006).

In the second step, all possible isotropic candidate structures described in the literature were examined; however, only two structures met the convergence criterion: the spatial power function (SP(POW)) and the MATERN function (SP(MATERN)).

The spatial power structure provides a direct generalization of the auto-regressive structure for equally spaced data and is assumed to be $f\,(d_{ij} ) = \rho ^{d_{ij}} $. The MATERN isotropic covariance function is given as $f\,(d_{ij} ) = \displaystyle{{(d_{ij} /2V)^V 2K_v (d_{ij} /\rho )} \over {\Gamma (V)}}$ where K _V is a modified Bessel function of the second kind of order V; Γ is the gamma function; ρ is a scalar parameter controlling the spatial range of correlation; and V is the ‘smoothness’ parameter, which allows greater flexibility for modelling the local spatial covariance. When V is small (V→0), the spatial process is assumed to be rough, and when V is large (V→∞), the process is assumed to be smooth (Littell et al. Reference Littell, Milliken, Stroup, Wolfinger and Schabenberger2006; Minasny & McBratney Reference Minasny and McBratney2007; SAS Institute Inc. 2008).

The NLMIXED procedure does not support the modelling of error structure directly; thus, the SAS macro %NLINMIX was used to refine the model to account for the chosen error structure. The restricted maximum likelihood was used along with the two expansion methods available in the SAS macro %NLINMIX (ZERO and EBLUP) to attempt to fit the model. The model adjusted in this way was compared with the fixed model and the random parameter model that was estimated using the previously described approach.

After selecting the model that best described the growth curves and included the previously estimated variance and covariance matrices, the hypotheses regarding one or more growth curve parameters could be tested. Thus, the difference between the growth curve parameters of the genotypes were tested, with the objective of simplifying the model (α=0·05).

Cross-validation

A cross-validation was performed to evaluate the goodness-of-fit of the chosen function. The original data set was randomly partitioned into two sub-data sets (training and test data), each containing 0·50 of the animals per genotype. The division of the original data set was conducted using PROC SURVEYSELECT from the SAS software (version 9.3). The first sub-data set (training data set) was used to fit the previously chosen best model (steps 1 and 2), and the second sub-data set (test data set) was used to evaluate the robustness of the model. Equation precision was measured with the coefficient of determination (R ²) between the observed and predicted values and the simultaneous F-test of the intercept and slope (intercept=0 and slope=1), whereas the accuracy of the equation was determined based on the concordance correlation coefficient (CCC) and the root-mean-square error prediction (RMSEP) and its decomposition into mean bias, systematic bias and random errors (Tedeschi Reference Tedeschi2006). The model evaluation was implemented using the Model Evaluation System v. 3.1.13 (http://nutritionmodels.tamu.edu/mes.html; verified 28 May 2013). Moreover, the parameters of the curve fit from the training data set were compared with the parameters of the curves fit from the original data set using 95% confidence intervals (IC_95%).

Impact of the number of weights on random parameters

One of the major advantages of the nonlinear mixed model methodology is the possibility of having individual growth parameters adjusted for each animal and testing to determine how some parameters differ from the population mean. Using the chosen model and the individual deviation of asymptotic and constant rate parameters, it was possible to evaluate the effect of the number of weights on the parameters estimated for each goat using descriptive graphic analysis.