Critical points on growth curves in autoregressive and mixed models

Adjusting autoregressive and mixed models to growth data fi ts discontinuous functions, which makes it diffi cult to determine critical points. In this study we propose a new approach to determine the critical stability point of cattle growth using a fi rst-order autoregressive model and a mixed model with random asymptote, using the deterministic portion of the models. Three functions were compared: logistic, Gompertz, and Richards. The Richards autoregressive model yielded the best fi t, but the critical growth values were adjusted very early, and for this purpose the Gompertz model was more appropriate.


Introduction
Growth functions have been used successfully to analyze responses of plants and animals to genetic selection and environmental variation.However, caution must be taken in adjusting these functions, since measurements in the same individual over time cannot be considered independent, since errors are serially correlated, and this information should be added to the model (Draper and Smith, 1998).Porter et al., 2010), Mendes et al. (2009), Mazzini et al. (2005), Bergamasco et al. (2001) and Calegario et al. (2005) have demonstrated the superiority of models with correlation structure in the errors compared to models that do not consider this correlation including.In all comparisons, autoregressive models were superior to models without autocorrelation.
The mixed effects model, in which a random effect parameter is added to the asymptote for each individual, accounts for increasing variances and serial correlations over time.Craig and Schinckel (2001) and Schinckel and Craig (2002) adjusted nonlinear mixed models to the weight of pigs, and concluded that this model outperformed the fi xed effects model.
Determining critical points in growth curves is important in agricultural and livestock work.For example, it allows one to adjust nutritional needs at the time when growth rate is at a maximum, and to determine the age at which an organism reaches maturity for vegetable or animal harvesting (Ersoy et al., 2006;Agudelo-Gómez et al., 2009 andTerra et al., 2010).
In autoregressive and mixed models the adjustment result is a discrete set of points, and it is impossible to determine critical points directly; it is necessary to use the deterministic part of the model, if it is well adjusted to the observed data.The aims of this study were to develop a methodology for determining critical points in autoregressive and mixed models logistic, Gompertz and Richards, and to compare the results with those obtained in models with independent errors.

Materials and Methods
The nonlinear model of fi xed effects that associate the y-response with x-time through the function F, without considering autocorrelation error, can be written as Model I: with θ = parameters vector, and  i = random error, independent, with normal distribution (0,  2


). Model I includes a deterministic part, F(x i , θ), and random errors  i considered to be independent, i.e. cov( i ,  j ) = 0 for i  j.Growth data obtained from measurements made on the same individual can, however, follow a model with dependent errors.If the i-th error  i is related to the previous one  i-1 , then replacing in the model (1)  i by   i-1 + Є i , where Є i is a random error, one obtains a fi rst-order autoregressive model, Model II: where ) and independent of  i , i. e. cov(Є i , Є j ) = 0 for i  j and cov( i , Є j ) = 0 ,  i, j.
Then, E(y From model I, y i = F(x i , θ) +  i , we can construct a mixed model that introduces a parameter  with random effects, allowing a variation of the upper asymptote , depending on the j-th individual, Model III: with  j = random effect of the j-th individual, with normal distribution (0,  2  ),  i,j = random error, independent, normally distributed (0,  2  ), independent of  j .The cov ( i,j ,  i',j' ) = 0 for i  i' or j  j', cov( i,j ,  i',j' ) =  2  for i = i' and j = j' and cov( i,j ,  j ) = 0  i, j; E(y i,j ) = F(x i , θ) e var(y i,j ) = [g(x i , θ)] 2  2  +  2  , which is separated into two components, variance between and within subjects.If g(x i , θ) is an increasing function with x then the model incorporates the structure of increasing variance with time.Furthermore, cov(y i,j , y i+k,,j ) = g(x i , θ) g(x i+k , θ)  2  , allowing the model to consider the correlation through time x.
With growth models, a variety of functions F(x i , ) can be employed, including logistic (Verhulst, 1845), Gompertz (Gompertz, 1825), Richards (Richards, 1959), von Bertalanffy (von Bertalanffy, 1957), and Brody (Brody, 1945).Only the Brody function has no infl ection point, and the Richards stands out for fl exibility in adjusting this point.The following functions will be considered: The residual autocorrelation may be verifi ed by statistical d Durbin-Watson, among other methods, according to Draper and Smith (1998), (9) where: ê i is the residual corresponding to the i-th observation, i = 1, ..., n = number of observations (x i ; y i ) in the model fi tting.The authors also present tables of critical values for probability levels.
The model adjustments can be checked with the usual criteria for comparing models, including RMS = residual mean square of the model; R 2 = square of the correlation coeffi cient between observed and estimated values by the model, as in Schinckel and Craig (2002); (10) where: SSR 1 = sum of squares of residuals of the model with smaller number of parameters, SSR 2 = sum of squares of residuals of the model with larger number of parameters, df 1 and df 2 = degrees of freedom associated with SSR 1 and SSR 2 , respectively; The mean prediction error in percentage is (11) where: y i = i-th observed value, ˆi y = i-th value estimated by the model.The MPE value allows one to check if the model under-or overestimates the observed data, depending on its sign (positive or negative, respectively); The Akaike Information Criterion (AIC), the corrected Akaike Information Criterion (AIC c ), and the Bayesian Information Criterion of Schwarz (BIC), respectively (12) (13) (14) where: SSR = Sum of squares of residuals of the model and p = number of model parameters.(Narinc et al., 2010).The information criteria evaluate whether the model adequately describes the studied population.The lower the value, the better the model.
The critical points of the functions, such as the infl ection point (pi) and deceleration asymptotic point (pda), are determined using the method developed in Mischan et al. (2011) for the model with independent errors, model I, for the logistic function.However, autoregressive models, model II, and mixed model, model III, lead to estimated values that are discrete variables defi ned only for the observed values of x i , i = 1, ..., n.Therefore we propose a new approach, which consists in the determination of the critical points through the function ( 15) the estimate of F(x i , θ) in ( 2) and ( 4), obtained by adjusting the models.
The abscissas of the points pi (x pi ) and pda (x pda ) are thus determined using the continuous function z, which should be well adjusted to the observed data.To verify the adequacy of the adjustment the following criteria are proposed: Residual mean square y i -z i 'nlmixed' procedure of SAS (Statistical Analysis System, versão 9.2) and the autoregressive model (model II) by the 'model' procedure using the macro '% ar'.The optimization algorithms used were Marquardt for the 'model' procedure and dual quasi-Newton for the 'nlmixed' procedure.Homoscedasticity of models was verifi ed by the Breusch-Pagan test (Breusch and Pagan, 1979) and the normality of the model errors was assessed by examination of residuals plotted against estimated values (Draper and Smith, 1998).
Growth data on weight of cows recorded monthly from birth (x = 0) to adult age, obtained in Piracicaba (22 o 42' S, 47 o 38' W, 546 m a.s.l.), state of São Paulo, Brazil, were used to illustrate the methodology.The observational data refer to fi ve Flemish breed (F) animals, nine from crossing Guersey × Gir (G) and 11 Holstein (H), evaluated for 60, 59, and 52 months respectively.

Results and Discussion
The following results compare models I, II, and III and the three growth functions logistic, Gompertz, and Richards fi tted to data of each animal (models I and II) and to animals grouped by race, or crossbreed (models I, II, and III).Considering the adjustment to each animal, the logistic and Gompertz functions fi t all 25 animals well, but convergence in the iterative method of fi tting the Richards function was obtained only for 21 animals in the case of model I and for 22 with model II.
Several authors have reported problems with the convergence process in the Richards adjustment.Brown et al. (1976), with weight-age data for female cattle, attributed the diffi culty of convergence in the Richards function mainly to the high negative value found for the residual correlations between the parameters  and , -0.97.In this study the average value for the residual correlations between the parameters was also -0.97 for the 21 animals in model I and -0.96 for the 22 in model II.This high correlation between the two parameters that account for the shape of the curve implies different combinations of the initial estimates in the iterative process of adjustment, which makes convergence diffi cult.
In the model fi tting to the animal group data, the adjustments obtained for the logistic function led to positive estimates a and c, and to negative estimates b, implying increasing curves with a positive superior asymptote and with an infl ection point.The same result was obtained in all adjustments of the Gompertz function, for which the estimates were a, b, and c > 0. For Richards, estimates a and c were always positive and b and d were of the same sign, characterizing increasing curves with a positive asymptote, with or without infl ection and minimum points.The parameter estimates are shown in Table 3 and the criteria used to compare models in Table 4.
With regard to the adjustment to the three groups of animals, the estimates of parameter  did not dif-( 17) where m = arbitrary number of points in the stabilization of growth.
The RMSz criterion verifi es the fi tting of the function throughout the growth interval of the individual, while the MPE m criterion verifi es this only in the fi nal stage.If deviations are equally distributed around z, then the MPE m value will be close to zero and the function will be well adjusted in this phase; this will give greater confi dence to the estimate of the critical point pda.The average deviation between the observed and estimated data in specifi c segments of the growth curve is one of the criteria suggested by Fitzhugh Jr. ( 1976) for model comparison.
The coordinates of the critical points pi (infl ection point) and pda (asymptotic deceleration point) were calculated as functions of the parameter estimates a, b, c, d.The Delta method of Oehlert (1992) was used to obtain the estimated variances of abscissas, (x pi ) and (x pda ); these are functions of the estimated parameters and the estimated variances and covariances of the estimated parameters: The formulae are in Tables 1 and 2.
If we consider increasing functions with a horizontal upper asymptote,  > 0, and  > 0, the Richards function, unlike the logistic and Gompertz, will have an infl ection point only for some combinations of signs of the parameters  and .Namely, if  < 0 and  < 0 the curve has an infl ection point; if  > 0 and  > 1 the curve has infl ection and minimum points.If  > 0 and 0 <  < 1 the curve has no infl ection point, only a minimum, whose coordinates are [ln (b) / c; 0].
The models with independent errors (model I) and with mixed effects (model III) were adjusted using the  Table 2 -Estimated variances for the abscissas of critical points of infl ection (x pi ), and of asymptotic deceleration (x pda ), for logistic, Gompertz, and Richards functions.
Sci. Agric.v.71, n.1, p.30-37, January/February 2014 fer greatly between models I and III, but are slightly larger in model II, autoregressive (Table 3).When the adjustment is made per animal, however, for logistic and Gompertz functions, in most cases, the estimates are not much different between the models I and II, but some cases are exception, with a reasonable decrease of up to 186 kg from the model I without autocorrelation to the autoregressive model.Mazzini et al. (2005) worked with Hereford males and found that the weighting of the data and the structure of autoregressive errors sub-stantially reduced the estimated values for asymptotes in the Gompertz and logistic functions, in the case of the adjustment to each animal.
Estimates c of the parameter  were very similar in all models and functions.As to the parameter , different models did not alter too much their estimated values, in the logistic and Gompertz functions, but in the case of Richards, the combination of estimated values of  and  led to functions with different shapes: the models I and III showed curves with similar aspect to the one  The approximate estimated standard errors of the estimates a and c were, on average, two to three times greater in the model II, autoregressive, in comparison with models I and III, for logistic, Gompertz and Richards functions.This increasing is expected, because the variance of the parameters is underestimated when we do not consider autoregressive structures (Souza, 1998).Concerning to the estimated standard error of the estimate b, the same can be said for the logistic and Gompertz functions, but when the Richards function is adjusted, the estimated standard errors of b and d are not comparable among different models, since the curves showed different growth patterns.
Durbin-Watson test for residuals (DW) in Table 4 was signifi cant (p < 0.05), showing positive correlation among residuals of fi tted models I and III.Model II, otherwise, was effi cient in the sense of incorporate the correlations to the predicted values, resulting in independent errors for the group F; for the groups G and H the values are very close to the superior critical signifi cance level.Model II was also effi cient in controlling heteroscedasticity, which was present in the model I. Furthermore, all values of r, the estimative of the autocorrelation parameter in model II, were signifi cant (p < 0.0001); thus one must preferably employ the autoregressive model in relation to models I and III.
All the goodness of fi t indicators, RMS, R 2 , AICc and BIC, showed that model II was the best one, followed by model III, and then by model I (Table 4).Mean square error in model II adjustment is, on average, four times smaller than the values in model I and 2.5 times smaller than in model III.F-statistic values in formula (10) were signifi cant (p < 0.05) and were performed to verify the contribution of Richards autoregressive model, which has one additional parameter, in relation to logistic and Gompertz models; this shows a signifi cant improvement in using this model.
Gompertz and Richards show better indicators of goodness of fi t than logistic (Table 4), i.e., smaller mean square errors (RMS) and Akaike's information criteria (AICc).R 2 values were similar in the adjustment of the three functions.On average, the logistic function presents lower estimates of asymptote, in relation to the other functions (Table 3).It seems to be a characteristic of logistics, underestimate mature weights and overestimate them in the early stages (Brown et al., 1976).This is probably due to the location of the point of infl ection, fi xed at one half of the asymptote.Figure 1 shows the fi t of Richards autoregressive and mixed models to the data of the Holstein race cows.
Second-order autoregressive models of logistic, Gompertz and Richards did not signifi cantly improve the adjustment, in relation to fi rst-order autoregressive model; residual mean squares were quite similar in both models (RMS-ar1/ RMS -ar2 = 1.01, in average).Mendes et al. (2009) worked with Hereford cows and observed a better fi t of the second-order autoregressive model Gompertz, in comparison with the fi rst-order, but the weighing was done from 15 to 15 days, whereas here it was monthly.
The adjustment of autoregressive Richards model to the group of animals, resulted in curves without infl ection points (Table 5).However, the adjustment of each animal led to some curves with infl ection and these are at the bottom of Table 5.Table 6 presents the MPE 20% and RMSz values.
In models I and III, the values of the abscissa of the infl ection point, on average of the three groups, decrease from the logistic function ( 16 months), where the ordinate is a/2, to Gompertz ( 11 months), where the ordinate is a/e (e = neperian number), and remain at this value in Richards, where the position of the point depends on the asymptote and the parameter  (Table 5).
In the autoregressive model II, the abscissas values are greater than in models I and III both in logistic function ( 20 months) as in Gompertz ( 13 months).Considering the group of 11 animals, where the Richards function was fi tted with infl ection point, it can be seen that the trend is the same of the 25 animals (lower part of Table 5); however in the autoregressive Richards model the abscissa is located very early (0.15 in average).The ratio ordinate of infl ection point to the asymptote, (y pi /a), that Fitzhugh Jr (1976)   rity at point of infl ection', was 0.372 in Richards models I and III, a similar value as the constant in Gompertz, 0.368, and smaller than 0.5 in the logistic function.In model II, however, the value is very early, the ratio y pi /a is only 0.080.This anticipation in the occurrence of the infl ection point of functions with fl exibility in its location is reported by Porter et al. (2010).
The asymptotic deceleration points, pda, were nearly the same for the three functions and models I and III ( 41 months, three races in average).A small increasing in this value was found for logistic and Gompertz in model II autoregressive.Concerning the analysis with 11 animals one sees a discrepant value for Richards autoregressive model (18.8 months, three races in average); at this age, the estimated weight is 278.1 in average, with ratio y pda /a equal to 0.529, a substantially lower value in comparison with the constants 0.908 in logistic and 0.847 in Gompertz functions.The weight that was reached in this point is only 53 % of the asymptote, which can not be indicated as a stability point, despite the better fi t of this function.As a comparison, a segmented regression (Portz et al., 2000) fi tted to data of the 25 animals led to an intersection point of two straight lines, one horizontal and other sloping, of 37.4 months.Figure 2 compares the position of pda points in Richards autoregressive and mixed models, for the fi t of the 11 animals' weights that had curves with infl ection points.
Almost all MPE 20% values are negative, meaning that z-function overestimates the observed values in the fi nal stage of growth (Table 6).This is more pronounced in model II in relation to models I and III, being the ratio M II / M III largest in logistics than in Gompertz.The situation is similar in the case of the adjustment to 11 animals, in addition that z-function is not so over-Sci.Agric.v.71, n.1, p.30-37, January/February 2014 estimated in the autoregressive Richards model.Using this criterion, MPE 20% , Richards function gives the best fi t, in the fi nal stages of growth, followed by Gompertz.The same conclusion may be attained when we used the RMSz criterion, that evaluates all the period of growth.Therefore, we can consider critical points determined by the autoregressive Richards model as well adjusted, and the relatively early values as a characteristic of the function, that is recognized as more fl exible with respect to the location of the infl ection point, and consequently of the asymptotic deceleration point.

Conclusions
With growth data of Flemish, Holstein and Guernsey-gir cows, the fi rst-order autoregressive model was effi cient in fi t the logistic, Gompertz and Richards functions, eliminating the residual autocorrelation.Mixed model, however, had a better fi t in relation to fi xed model with independent errors, but was not effi cient in correcting the residual autocorrelation effect.Concerning the function, Richards fi ts better than logistic and Gompertz.In the determination of critical points in growth curves adjusted with models with a fi xed and a random part, as the autoregressive and mixed models here considered, the deterministic function must be employed.The criteria proposed to verify the fi t of the deterministic function, i.e., residual mean square for this function, and mean prediction error in the fi nal observations, have allowed the comparison of the This showed that Richards fi ts better than Gompertz, and this is better than the logistic.
The values obtained to estimate the abscissa of the deceleration asymptotic point in logistic and Gompertz autoregressive model can be used as a stability point, but the value obtained with Richards autoregressive model was a very early stability point.
number of observed data and p = number of parameters of the model; MPE criterion defi ned in (11), adapted to the function z in the stabilization stage of growth, Sci.Agric.v.71, n.1, p.30-37, January/February 2014 no explicit formula for x pda and y pda for the Richards function.

Figure 1 -
Figure 1 -First-order autoregressive (M II) and mixed (M III) Richards models fi tted to 11 Holstein (H) cow weights.yo = observed values, ye = estimated values, z = deterministic part of the model.

Figure 2 -
Figure 2 -Asymptotic deceleration points, pda, fi tted to data of three groups of cows, Flemish (F), Guernsey × Gir (G) and Holstein (H), through autoregressive (M II, y pda /a = 0.529) and mixed (M III, y pda /a = 0.848) Richards models.11 animals in each group.

Table 1 -
Coordinates of the critical points of infl ection (x pi ; y pi ) and of asymptotic deceleration (x pda ; y pda ) for logistic, Gompertz, and Richards functions.

Table 3 -
presented by the Gompertz and logistic functions, with negative estimates b and d, and the model II has presented positive estimates b and d, but with d < 1, which determines a curve without infl ection point.

Table 5 -
Estimates of infl ection (pi) and asymptotic deceleration (pda) points for logistic, Gompertz, and Richards functions fi tted by models I (with independent errors), II (autoregressive) and III (mixed), and estimated standard errors (ep) for the abscissas of the points.Flemish (F), Guernsey × Gir (G), and Holstein (H) cow groups.

Table 6 -
Comparison criteria for z-function fi tting of logistic, Gompertz, and Richards functions with models I (with independent errors), II (autoregressive), and III (mixed): mean prediction error for the latest 20 % observations (MPE 20% ) and residual mean square of z-function (RMSz).Flemish (F), Guernsey × Gir (G), and Holstein (H) cow groups.