SciELO - Scientific Electronic Library Online

 
vol.71 issue1Forage yield and nutritive value of Panicum maximum genotypes in the Brazilian savannahMorphological phenotypic dispersion of garlic cultivars by cluster analysis and multidimensional scaling author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

Share


Scientia Agricola

On-line version ISSN 1678-992X

Sci. agric. (Piracicaba, Braz.) vol.71 no.1 Piracicaba Jan./Feb. 2014

http://dx.doi.org/10.1590/S0103-90162014000100004 

BIOMETRY, MODELING AND STATISTICS

 

Critical points on growth curves in autoregressive and mixed models

 

 

Sheila Zambello de Pinho*; Lídia Raquel de Carvalho; Martha Maria Mischan; José Raimundo de Souza Passos

São Paulo State University/IBB - Dept. of Biostatistics, C.P. 510 - 18618-970 - Botucatu, SP - Brazil

 

 


ABSTRACT

Adjusting autoregressive and mixed models to growth data fits discontinuous functions, which makes it difficult to determine critical points. In this study we propose a new approach to determine the critical stability point of cattle growth using a first-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 fit, 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 fixed 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 and Terra 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 fixed 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(xi, θ), 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 first-order autoregressive model, Model II:

where

θ1 = parameters vector = [θ' ρ]', ρ = correlation measurement between errors, -1 < ρ < 1, and Єi = random error, independent, normally distributed (0, σ2Є) and independent of εi, i. e. cov(Єi, Єj) = 0 for i ≠ j and cov(εi, Єj) = 0 , ∀ i, j.

Then, E(yi) = F(xi, θ), var(yi) = ρ2 σ2ε + σ2Є and cov(yi, yj) = ρ2 cov(εi-1, εj-1) for i ≠ j.

From model I, yi = F(xi, θ) + ε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(yi,j) = F(xi, θ) e var(yi,j) = [g(xi, θ)]2 σ2Δ + σ2ξ, which is separated into two components, variance between and within subjects. If g(xi, θ) is an increasing function with x then the model incorporates the structure of increasing variance with time. Furthermore, cov(yi,j, yi+k,,j) = g(xi, θ) g(xi+k, θ) σ2Δ, allowing the model to consider the correlation through time x.

With growth models, a variety of functions F(xi, θ) 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 inflection point, and the Richards stands out for flexibility in adjusting this point. The following functions will be considered:

(a) Logistic

(b) Gompertz

(c) Richards

δ = 1/(1- υ), θ = [α β γ υ]', or θ = [α β γ δ]', with α > 0, β < 0, γ > 0 and υ > 1 (δ < 0), or α > 0, β > 0, γ > 0 and υ < 1 (δ > 0).

The residual autocorrelation may be verified by statistical d Durbin-Watson, among other methods, according to Draper and Smith (1998),

where: êi is the residual corresponding to the i-th observation, i = 1, ..., n = number of observations (xi; yi) in the model fitting. 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; R2 = square of the correlation coefficient between observed and estimated values by the model, as in Schinckel and Craig (2002);

where: SSR1 = sum of squares of residuals of the model with smaller number of parameters, SSR2 = sum of squares of residuals of the model with larger number of parameters, df1 and df2 = degrees of freedom associated with SSR1 and SSR2, respectively;

The mean prediction error in percentage is

where: yi = i-th observed value, ŷi = 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 (AICc), and the Bayesian Information Criterion of Schwarz (BIC), respectively

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 inflection 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 defined only for the observed values of xi, i = 1, ..., n. Therefore we propose a new approach, which consists in the determination of the critical points through the function

the estimate of F(xi) in (2) and (4), obtained by adjusting the models.

The abscissas of the points pi (xpi) and pda (xpda) 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 yi - zi

with n = number of observed data and p = number of parameters of the model; MPE criterion defined in (11), adapted to the function z in the stabilization stage of growth,

where m = arbitrary number of points in the stabilization of growth.

The RMSz criterion verifies the fitting of the function zi = F(xi,) throughout the growth interval of the individual, while the MPEm criterion verifies this only in the final stage. If deviations are equally distributed around z, then the MPEm value will be close to zero and the function will be well adjusted in this phase; this will give greater confidence to the estimate of the critical point pda. The average deviation between the observed and estimated data in specific 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 (inflection 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, (xpi) and (xpda); these are functions of the estimated parameters and the estimated variances and covariances of the estimated parameters: fbb = (b), fcc = (c), fdd = (d), fbc = côv (b,c), fbd = côv (b,d), fcd = côv (c,d). 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 inflection point only for some combinations of signs of the parameters β and δ. Namely, if β < 0 and δ < 0 the curve has an inflection point; if β > 0 and δ > 1 the curve has inflection and minimum points. If β > 0 and 0 < δ < 1 the curve has no inflection 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 '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 verified 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º42' S, 47º38' W, 546 m a.s.l.), state of São Paulo, Brazil, were used to illustrate the methodology. The observational data refer to five 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 fitted 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 fit all 25 animals well, but convergence in the iterative method of fitting 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 difficulty 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 difficult.

In the model fitting 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 inflection 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 inflection 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 differ 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 substantially 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 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 inflection point.

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 significant (p < 0.05), showing positive correlation among residuals of fitted models I and III. Model II, otherwise, was efficient 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 significance level. Model II was also efficient 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 significant (p < 0.0001); thus one must preferably employ the autoregressive model in relation to models I and III.

All the goodness of fit indicators, RMS, R2, 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 significant (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 significant improvement in using this model.

Gompertz and Richards show better indicators of goodness of fit than logistic (Table 4), i.e., smaller mean square errors (RMS) and Akaike's information criteria (AICc). R2 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 inflection, fixed at one half of the asymptote. Figure 1 shows the fit 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 significantly improve the adjustment, in relation to first-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 fit of the second-order autoregressive model Gompertz, in comparison with the first-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 inflection points (Table 5). However, the adjustment of each animal led to some curves with inflection and these are at the bottom of Table 5. Table 6 presents the MPE20% and RMSz values.

 

 

In models I and III, the values of the abscissa of the inflection 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 fitted with inflection 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 inflection point to the asymptote, (ypi/a), that Fitzhugh Jr (1976) called the 'degree of maturity at point of inflection', 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 ypi/a isonly 0.080. This anticipation in the occurrence of the inflection point of functions with flexibility 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 ypda/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 fit of this function. As a comparison, a segmented regression (Portz et al., 2000) fitted 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 fit of the 11 animals' weights that had curves with inflection points.

 

 

Almost all MPE20% values are negative, meaning that z-function overestimates the observed values in the final 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 overestimated in the autoregressive Richards model. Using this criterion, MPE20%, Richards function gives the best fit, in the final 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 flexible with respect to the location of the inflection point, and consequently of the asymptotic deceleration point.

 

Conclusions

With growth data of Flemish, Holstein and Guernsey-gir cows, the first-order autoregressive model was efficient in fit the logistic, Gompertz and Richards functions, eliminating the residual autocorrelation. Mixed model, however, had a better fit in relation to fixed model with independent errors, but was not efficient in correcting the residual autocorrelation effect. Concerning the function, Richards fits better than logistic and Gompertz.

In the determination of critical points in growth curves adjusted with models with a fixed and a random part, as the autoregressive and mixed models here considered, the deterministic function must be employed. The criteria proposed to verify the fit of the deterministic function, i.e., residual mean square for this function, and mean prediction error in the final observations, have allowed the comparison of the functions. This showed that Richards fits 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.

 

References

Agudelo-Gómez, D.; Hurtado-Lugo, N.; Cerón-Muñoz, M.F. 2009. Growth curves and genetic parameters in Colombian buffaloes (Bubalus bubalis Artiodactyla, Bovidae). Revista Colombiana de Ciencias Pecuarias 22:178-188.         [ Links ]

Bergamasco, A.F.; Aquino, L.H.de; Muniz, J.A. 2001. Non linear models fitting to growth data in females of the Holstein heifers. Ciência e Agrotecnologia 25:235-241 (in Portuguese, with abstract in English).         [ Links ]

Bertalanffy, L. von. 1957. Quantitative laws for metabolism and growth. The Quarterly Review of Biology 32:217-231.         [ Links ]

Breusch, T.S.; Pagan, A.R. 1979. A simple test for heteroscedasticity and random coefficient variation. Econometrica 47:1287-1294.         [ Links ]

Brody, S. 1945. Bioenergetics and Growth. Reinhold, New York, NY, USA.         [ Links ]

Brown, J.E.; Fitzhugh Jr, H.A.; Cartwright, T.C. 1976. A comparison of nonlinear models for describing weight-age relationships in cattle. Journal of Animal Science 42:810-818.         [ Links ]

Calegario, N.; Calegario, C.L.L.; Maestri, R.; Daniels, R. 2005. Improving fitting quality of forest biometric models by applying the generalized nonlinear model theory. Scientia Florestalis 69:38-50 (in Portuguese, with abstract in English).         [ Links ]

Craig, B.A.; Schinckel, A.P. 2001. Nonlinear mixed effects model for swine growth. The Professional Animal Scientist 17:256-260.         [ Links ]

Draper, N.R.; Smith, H. 1998. Applied Regression Analysis. 3ed. John Wiley, New York, NY, USA. 706p.         [ Links ]

Ersoy, I.E.; Mendes, M.; Aktan, S. 2006. Growth curve establishment for American Bronze turkeys. Archiv Tierzucht 49:293-299.         [ Links ]

Fitzhugh Jr, H.A. 1976. Analysis of growth curves and strategies for altering their shape. Journal of Animal Science 42:1036-1051.         [ Links ]

Gompertz, B. 1825. On the nature of the function expressive of the law of human mortality, and on a new method of determining the value of life contingencies. Philosophical Transactions of the Royal Society 182:513-585.         [ Links ]

Mazzini, A.R.A; Muniz, J.A.; Silva, F.F.; Aquino, L.H. 2005. Growth curve for hereford males: heterocedasticity and autoregressives residuals. Ciência Rural 35:422-427 (in Portuguese, with abstract in English).         [ Links ]

Mendes, P.N.; Muniz, J.A.; Silva, F.F.; Mazzini, A.R.A.; Silva, N.A.M. 2009. Analysis of the difasics growth curve of Hereford females by the Gompertz non-linear function. Ciência Animal Brasileira 10:454-461 (in Portuguese, with abstract in English).         [ Links ]

Mischan, M.M.; Pinho, S.Z.; Carvalho, L.R. 2011. Determination of a point sufficiently close to the asymptote in nonlinear growth functions. Scientia Agricola 68:109-114.         [ Links ]

Narinc, D.; Karaman, E.; Firat, M.Z.; Aksoy, T. 2010. Comparison of non-linear growth models to describe the growth in Japanese Quail. Journal of Animal and Veterinary Advances 9:1961-1966.         [ Links ]

Oehlert, G.W. 1992. A note on the delta method. The American Statistician 46:27-29.         [ Links ]

Porter, T.; Kebreab, E.; Darmani Kuhi, H.; Lopez, S.; Strathe, A.B.; France, J. 2010. Flexible alternatives to the Gompertz equation for describing growth with age in turkey hens. Poultry Science 89:371-378.         [ Links ]

Portz, L.; Dias, C.T.S.; Cyrino, J.E.P. 2000. A broken-line model to fit fish nutrition requirements. Scientia Agricola 57:601-607.         [ Links ]

Richards, F.J. 1959. A flexible growth function for empirical use. Journal of Experimental Botany 10:290-300.         [ Links ]

Schinckel, A.P.; Craig, B.A. 2002. Evaluation of alternative nonlinear mixed effects models of swine growth. The Professional Animal Scientist 18:219-226.         [ Links ]

Souza, G.S. 1998. Introduction to Models of Linear and Nonlinear Regression = Introdução aos Modelos de Regressão Linear e Não Linear. Embrapa-SPI/Embrapa-SEA, Brasília, DF, Brasil. 489 p.         [ Links ]

Terra, M.F.; Muniz, J.A.; Savian, T.V. 2010. Fitting the logistic and Gompertz models to growth data of fruits of date palm-dwarf (Phoenix roebelenii O'BRIEN) = Ajuste dos modelos logístico e Gompertz aos dados de crescimento de frutos de tamareira-anã (Phoenix roebelenii O'BRIEN). Magistra 22:01-07.         [ Links ]

Verhulst, P.-F. 1845. Mathematical researches into the law of population growth increase. Nouveaux Mémoires de l'Academie Royale des Sciences et Belles-Lettres de Bruxelles 18:1-42.         [ Links ]

 

 

Received April 15, 2013
Accepted July 15, 2013

 

 

Edited by: Thomas Kumke
* Corresponding author <sheila@ibb.unesp.br>

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License