Bayesian Analysis for the Comparison of Nonlinear Regression Model Parameters: an Application to the Growth of Japanese Quail

ABSTRACT

This paper discusses the Bayesian approach as an alternative to the classical analysis of nonlinear models for growth curve data in Japanese quail. A Bayesian nonlinear modeling method is introduced and compared with the classical nonlinear least squares (NLS) method using three non-linear models that are widely used in modeling the growth data of poultry. The Gompertz, Richards and Logistic models were fitted to 499 Japanese quail weekly averaged body weight data. Normal prior was assumed for all growth curve parameters of the models with assuming Jeffreys' non-informative prior for residual variances. Models were compared based on the Bayesian measure of fit, deviance information criterion (DIC), and our results indicated the better fit of Gompertz and Richards models than the Logistic model to our data. Moreover, the parameter estimates of the models fitted by both approaches showed only small differences.

Keywords:
Gompertz; Logistic; Richards; non-linear; Bayesian

INTRODUCTION

While a large number of linear functions are adequately used to model a wide variety of relationships between variables, many biological phenomena require non-linear functions, in which the response varies as a function of time. Time-related changes of a phenomenon are of particular importance in a wide range of disciplines such as biology, agriculture, economics, medicine, crop science, etc. Growth curves illustrating these changes allow the data to be summarized by a few number of parameters known as growth curve parameters.

There has been a great deal of interest in modeling the growth of poultry. Most studies have successfully fitted the Gompertz, Richards, and Logistic functions to Japanese quail data (Marks, 1978; Anthony et al., 1991; Narinc et al., 2009, 2010). On the other hand, several other studies have modeled the growth of other species of poultry (Tzeng & Becker, 1981; Emmans, 1995;Akbas & Oguz, 1998; Aggrey, 2002;Ahmadi & Golian, 2008).

Growth curves are generally fitted by nonlinear regression or linear regression if the model can be linearized by transformation. However, a linear form of most of the widely used growth models does not exist (Blasco et al., 2003). Nonlinear functions are particularly suitable for modeling growth data, since predictions outside the range of data set can be obtained more reliably than by linear models, and few parameters having a biological interpretation can be used to describe the entire growth process (Vuori et al., 2006). Nonlinear models are more complicated to solve than linear models, and therefore, many different algorithms have been developed to obtain the estimates of model parameters. More recently, an alternative to the classical approach, the Bayesian approach, has received much attention as a tool for estimating the parameters of a variety of nonlinear models, including complex growth curve models (Zhang et al., 2007).

In the Bayesian approach, as an alternative to the frequentist or classical approach, the assumption of normality of the data set is not a necessary condition, and the inferences on the model parameters are based on their posterior distribution. In general, they produce more accurate estimates, with greater representativeness of the statistical models and have been shown to be easy to apply and understand. It is a tool with great potential, as it can consider the uncertainty of all the parameters in a model. Bayesian methods also provide a natural and principled way to incorporate prior information about unknown parameters with the data in hand, which can improve the accuracy of the results or predictions. The key point of the Bayesian approach is that it treats parameters as random variables, whereas they are fixed constants in the classical approach. Therefore, Bayesian approach may be an alternative to the frequentist methodology for fitting nonlinear models to the growth data of Japanese quail. Although the Bayesian approach is powerful and can be used in a variety of analytic models, the difficulties of programming and computational demands have discouraged many researchers from using it (Zhang et al., 2007). However, the use of the Markov Chain Monte Carlo (MCMC) methods have allowed the researchers to numerically obtain the solution of complex posterior density integration and made the Bayesian methods popular among the researchers from different disciplines.

Although Bayesian methods have become well established in the animal science literature and have been successfully applied to linear models (Theobald et al., 1997; Firat et al., 1997a; Firat et al., 1997b; Firat, 2001), to our knowledge, none of the studies employed Bayesian methods to estimate the growth curve parameters of Japanese quail in nonlinear models. This study aimed at fitting three commonly used nonlinear models (Gompertz, Richards and Logistic models) to describe the growth curves of Japanese quail. Growth curve parameters were estimated using the classical (NLS) approach, as well as the Bayesian method, allowing estimating the joint posterior distribution of growth curve parameters. Models were compared using the deviance information criterion (DIC) that permits the comparison of the overall predictive ability.

MATERIAL AND METHODS

Animals

The data set, which was first published in Narinc et al. (2010), were formed the basis of this study. The study was conducted at the Poultry Breeding Unit, Animal Science Department, Faculty of Agriculture, Akdeniz University, Turkey. Averages of weekly recorded body weight measurements belonging to 499 Japanese quail were used in model fitting.

Models' Formulation

A growth model for a single experimental unit was assumed as:

v

where is the observed weight, n is total number of observations; is a vector of unknown parameters, t_j is the time of which the j ^th observation was taken, is independent random error of y_j , and is one of the different types of growth curves. For the analysis of our data set, the following three nonlinear functions, which are frequently used for the description of growth curves of Japanese quail, were considered:

Gompertz model:

where

Richards model:

where

Logistic model:

where

In the above models, 0 stands for the asymptotic weight, 1 is a constant without biological meaning, 2 is the maturity index as an expression of the growth rate relative to mature weight, and 3 is the shape parameter defining the position of the point of inflection (Aggrey, 2002; Kizilkaya et al., 2006).

Then, were substituted in the growth model and denoted y= and ,their likelihood functions were obtained, as follows:

where ,

Bayesian Inference

The following distributions of the three growth model priors ar proposed:

Gompertz model:

Then .

Hence, the prior distribution is

Richards model:

, and , and b3 are independent, the vector .

Hence, the prior distribution is:

Logistic model:

where are independent, the vector .

Hence, the prior distribution is:

The prior distribution for the error variances are assigned a scaled inverted chi-square distribution with hyper-parameters v and s2

where u is the degree of belief and s2 is the scaling factor of the inverted chi-square distribution.

In practice, the specification of hyper-parameters may be difficult, and therefore we can choose the values of hyper parameters to obtain non-informative priors. For example, setting u and s2 to zero will yield Jeffreys' prior for s2 and thus we will have density .

The joint posterior distributions of all the unknown parameters for a single experimental unit i are as follows:

Gompertz model:

Richards model:

Logistic model:

The joint posterior distributions described in (1)-(3) are complex, and the implementation of the Bayesian procedure requires MCMC sampling techniques, and in particular, the Gibbs sampling and Metropolis-Hastings (MH) algorithm. The Gibbs sampler generates posterior samples using the conditional densities of the model parameters. The full conditional posterior distributions of each model parameter and error variance for running the Gibbs sampler and Metropolis-Hastings are needed to derive. These full conditional distributions can be obtained from algebraic and matrix manipulations of the joint posterior distributions in equations (1)-(3). The full conditional distribution of error variance, was proportional to inverted distribution,

Suppose we take the Richards model as an example. The full conditionals for the unknown parameters from (2) are as follows:

Conditional posterior distributions of the growth curve parameters (b_0, b_1, b_2, b₃) of models in equations (5)-(8) do not have a closed form. Therefore, these parameters require implementation of a Metropolis algorithm, which is used to sample from the appropriate conditional posterior densities given in (5)-(8). However, the error variance, s², has a scaled inverse Chi-square conditional posterior distribution in equation (4) and Gibbs sampling can be used to draw samples from this recognizable conditional distribution. Metropolis steps can be embedded in the Gibbs scheme in order to draw samples from the non-standard statistical distributions.

Goodness of Fit

There exist a variety of methodologies to compare several competing models for a given data set and to select the one that best fits the data. In this study, models were compared using the Bayesian measure of fit, the deviance information criterion (DIC), which is a tool that is used for model assessment and provides a Bayesian alternative to Akaike's information criterion (AIC) and the Bayesian information criterion (BIC). This statistic takes into account the number of unknown parameters in the model and is intended as a generalization of the AIC. Let be the parameters of the i'th model, the DIC is computed as follows (Forni et al., 2009):

where,

and

In model comparison, the smaller the DIC model fits better to the data set (Spiegelhalter et al., 2002; Forni et al., 2009). In practical applications of DIC, Spiegelhalter et al. (2003) stated that just reporting the model with the lowest DIC could be misleading if the difference in DIC is small, for example less than 5, and the models make very different inferences. DIC gives a clear conclusion to support the null hypothesis or the alternative hypothesis similar to the Bayes factor, BIC, and AIC.

Statistical Analysis

Parameter estimates of the three non-linear models were obtained by two methods:

i) Frequentist or classical: A non-linear least square (NLS) problem is an unconstrained minimization problem, where the objective function is non-linear in terms of model parameters. Unlike the linear least square estimation, closed form solutions to the k normal equations cannot be obtained, and therefore, iterative methods such as Gauss-Newton are required to minimize the sum of squares. Therefore, for estimating the parameters of growth functions, the Gauss-Newton algorithm was utilized and the NLIN procedure of SAS 9.2 software (SAS Ins., 2009) was used.

ii) Bayesian: Posterior samples of the model parameters were obtained by the WinBUGS program, 1.4.2 version (Lunn et al., 2000) through the R2WinBUGS package of R program (R Development Core Team, 2010). A normal distribution with mean 0 and variance 10,000, N(0, 10,000), was used as the prior distribution for each of the growth curve parameters in models (), and a gamma distribution Ga(0.001,0.001) was used as prior for the precision parameters. Both N(0, 10,000) and Ga(0.001, 0.001) are considered "vague". The distribution Ga(0.001, 0.001) (a proper distribution) is very close to the commonly used Jeffreys' prior for the precision parameter, an improper distribution that does not integrate to one. Using this "barely" proper distribution, it is possible to avoid improper posteriors under the Jeffreys' prior. The distribution N(0, 10,000), with a mean 0 and standard deviation of 100, is practically flat. Vague prior distributions were used because we had no information about these model parameters. In all models, chain lengths of 500,000 cycles were considered with a burn-in period of first 50,000 cycles. Every 225th sample is retained in the next 450,000 samples to reduce the auto-correlation and in total 2,000 samples were used to estimate the features of marginal distributions.

RESULTS AND DISCUSSION

Parameter estimates, their approximate standard errors and 95% confidence intervals obtained from the classical approach (NLS) for the three growth functions are given in Table 1 for completeness. Also, Table 2 provides the summary statistics of Bayesian analysis for the fitted functions.

Table 1
NLS estimation results for the three growth functions

		Estimate	Approximate Std. Error (SE)	95% Approximate Confidence Interval (CI)
				2.50%	97.50%
Gompertz	b0	222.10	1.584	217.700	226.500
	b1	3.312	0.044	3.190	3.435
	b2	0.080	0.001	0.077	0.084
	AIC	-2.010
	BIC	-2.172
Richards	b0	222.00	3.191	211.800	232.100
	b1	0.010	0.237	-0.743	0.763
	b2	0.080	0.004	0.066	0.094
	b3	0.003	0.071	-0.223	0.229
	AIC	-0.014
	BIC	-0.231
Logistic	b0	201.90	5.982	185.300	218.600
	b1	12.824	2.098	6.999	18.649
	b2	0.139	0.011	0.108	0.169
	AIC	24.378
	BIC	24.215

• AIC: Akaike's information criterion; BIC: Bayesian information criterion

Table 2
Posterior summary statistics for the three growth functions

		Median	Mean	SD	MCError	95% Credible Interval (BCI)
						2.50%	97.50%
Gompertz	b0	222.10	222.14	2.14	0.080	217.80	226.60
	b1	3.312	3.312	0.057	0.002	3.195	3.429
	b2	0.080	0.080	0.002	0.000	0.077	0.083
	s2	0.653	1.031	1.871	0.044	0.198	4.059
	DIC	22.3
Richards	b0	221.80	221.86	2.19	0.071	217.70	226.31
	b1	0.014	0.015	0.003	0.000	0.012	0.021
	b2	0.080	0.080	0.002	0.000	0.077	0.084
	b3	0.004	0.005	0.001	0.000	0.004	0.007
	s2	0.646	1.077	1.772	0.047	0.200	4.766
	DIC	22.3
Logistic	b0	202.00	202.58	7.26	0.254	189.70	217.30
	b1	12.630	12.586	1.341	0.186	9.491	14.940
	b2	0.138	0.138	0.009	0.001	0.117	0.156
	s2	23.870	34.076	34.914	1.238	8.060	117.805
	DIC	46.9

• SD: Standard deviation; MCError: Monte Carlo Error; DIC: Deviance information criterion

Gompertz function has the best fit to the data in terms of model selection criteria, AIC and BIC, followed by Richards and Logistic functions using classical approach, NLS. On the other hand, Figure 1 depicts the fit of three growth functions to the actual values within Bayesian context. From the figure, it can clearly be seen that all models are accurate in prediction of the actual body weights. The DIC statistics of the aforementioned models indicated that both the Gompertz and Richards functions fits the data exactly the same and also are better than the logistic function (Table 2). Comparing the goodness of fit of various growth functions to quail data, numerous studies have shown that Gompertz and Richards functions reflect the age-weight relationship of Japanese quails well (Tzeng & Becker, 1981; Anthony et al., 1991; Akbas & Oguz, 1998; Mignon-Grasteau et al., 1999; Sezer &Tarhan, 2005; Alkan et al., 2009; Narinc et al., 2010).

Figure 1
Actual and Fitted body weights (g) for Gompertz, Richards and Logistic functions via Bayesian approach.

Similar parameter estimates were obtained when classical and Bayesian approaches were considered (Tables 1 and 2). The standard deviations (SD) of the posterior distribution are analogous to the approximate standard error (SE) estimates of the parameters which can be produced by default in any statistical package. In this respect, SDs (Table 2) were used as equivalent estimates of the SEs (Table1). For the Gompertz function, the Bayesian analysis resulted higher SDs of the parameter estimates than their corresponding SEs. The Bayesian estimation of the Richards function produced smaller SDs for all the parameters. While the posterior SD of b₀ is higher than its SE in NLS estimation, smaller SEs are obtained for b₁ and b₂ when fitting the logistic function. In the Bayesian context, the classical confidence intervals (CI) are replaced by the Bayesian credible intervals (BCI). Both CI and BCI for the Gompertz function parameters were similar, while BCIs are narrower than CIs for Logistic and Richards function parameters (Tables 1 and 2). On the other hand, CIs of both b₁ and b₃ parameters of Richards function include zero, and are considered to be statistically insignificant. Non-significance of any of the estimated parameters of a growth function might have several reasons (Fekedulegn et al., 1999): 1) one or more parameters in the function may not be useful, or model can be reparameterized to involve less parameters; 2) the growth data are not adequate to estimate all parameters in the function; or 3) the assumptions of the function do not conform with the biological system being modeled. Moreover, it was reported that convergence problems may occur in fitting Richards function (Darmani Kuhi et al., 2003; Karkach, 2006). Although Logistic function led to statistically significant parameter estimates using NLS approach, as can be seen by a better fit of Richards function than the Logistic function, significant parameter estimates do not always guarantee a good fit to the data, and therefore this should not be the only criterion to select the best fitting model for the data in hand (Darmani Kuhi et al., 2003). Contrary to NLS method, BCIs of the parameters in Richards function do not include zero. This is due to the power of Bayesian approach in estimating the parameters of complex functions, particularly with small data sets.

Although the parameters of the functions are represented with the same letters, they do not have the same biological meaning, except the asymptotic weight parameter, b₀, and comparison of the parameter estimates across the models is not reasonable. The asymptotic weight parameter estimates obtained from fitting Gompertz, Richards and Logistic functions are 222.10, 221.80 and 202.00, respectively. These values of b₀ parameter are higher than Anthony et al. (1991) have reported. Similar results have been reported by Akbas & Oguz (1998), Kizilkaya et al. (2005) and Alkan et al. (2009). The differences in the parameter estimates of growth functions may occur among different studies as a result of having different genotype and environmental conditions.

Gompertz and Logistic functions have a priori defined growth shapes with inflection point at 37% and 50% of the asymptote, respectively (Grimm & Ram, 2009). In the Richards function, an additional parameter, b₃, allows the flexibility in the inflection point. Moreover, the Gompertz and Logistic growth functions can be derived from the Richards function. When the b₃ parameter in Richards function approaches 0, Gompertz function is approximated, where as if b₃ parameter approaches 1, the Logistic function is obtained (Sezer & Tarhan, 2005;Kizilkaya et al., 2006). When choosing an appropriate function to model growth in biological systems, researchers need to pick a model that is flexible with the least complexity among the available functions. For instance, the Gompertz and Logistic functions are simple and were shown to fit well to short time series such as growth data of animal species (Darmani Kuhi et al., 2010). Compared to Gompertz and Logistic functions, Richards function is more complex, has an additional parameter and can fit well to the complex patterns but it is difficult to fit and requires a long time series (Karkach, 2006).

CONCLUSIONS

In this paper, we considered the estimation of the growth curve parameters for nonlinear models using classical and Bayesian methods. Both classical and Bayesian methods of estimation are well established and it is not necessary to justify why one or the other is preferred (Blasco, 2001). The Bayesian methodology can be applied without restriction to the growth functions used in this study. However, Bayesian analyses are sensitive to the choice of priors and different priors may result in different posterior expectations of the parameters. Since the parameters of the normal priors (hyper parameters) were chosen so as to yield nearly a flat distribution, classical and Bayesian estimates did not differ much for the nonlinear curve parameters. Further study may be carried out to examine the effects of different prior specifications in nonlinear functions and to estimate their parameters in Bayesian context. However, the Bayesian method consistently predicts narrower credible intervals and gives more accurate results than those from the classical approach. When using the Bayesian methodology, it is possible to proceed to the comparisons among these parameters without having to resort to asymptotic procedures and incoherent results through frequentist theories. It is also important to note that use of these procedures need not be confined to any specific class of models. Rather, they can be applied to a wide range of model determination problems.

ACKNOWLEDGEMENT

We would like to express our gratitude to Dr. Tulin Aksoy (Akdeniz University, Faculty of Agriculture, Department of Animal Science, Antalya-TURKEY) for supplying the data, which was collected as part of a research project supported by Akdeniz University Scientific Research Projects Management Unit (project number: 2005.02.0121.005).

REFERENCES

• Aggrey SE. Comparison of three nonlinear and spline regression models for describing chicken growth curves. Poultry Science 2002;81:1782-1788.
• Ahmadi H, Golian A. Non-linear hyperbolastic growth models for describing growth curve in classical strain of broiler chicken. Research Journal of Biological Sciences 2008;3(11):1300-1304.
• Akbas Y, Oguz I. Growth curve parameters of lines of Japanese quail (Coturnix coturnix japonica), unselected and selected for four-week body weight. Archiv Fur Geflugelkunde 1998;62(3):104-109.
• Alkan S, Mendes M, Karabag K, Balcioglu MS. Effect of short-term divergent selection for 5-week body weight on growth characteristics of Japanese quail. Archiv Fur Geflugelkunde 2009;73(2):124-131.
• Anthony NB, Emmerson DA, Nestor KE, Bacon WL. Comparison of growth curves of weight selected populations of turkeys, quail and chickens. Poultry Science 1991;70:13-19.
• Blasco A, Piles M, Varona L. A Bayesian analysis of the effect of selection for growth rate on growth curves in rabbits. Genetics Selection Evolution 2003;35:21-41.
• Blasco A. The Bayesian controversy in animal breeding. Journal of Animal Science 2001;79:2023-2046.
• Darmani Kuhi H, Kebreab E, Lopez S, France J. An evaluation for different growth functions for describing the profile of live weight with time (age) in meat and egg strains of chicken. Poultry Science 2003;82:1536-1543.
• Darmani Kuhi H, Porter T, Lopez S, Kebreab E, Strathe AB, Dumas A, et al. A review of mathematical functions for the analysis of growth in poultry. World's Poultry Science Journal 2010;66:227-239.
• Emmans GC. Problems in modeling the growth of poultry. World's Poultry Science Journal 1995;51:77-89.
• Fekedulegn D, Mac Siurtain MP, Colbert JJ. Parameter Estimation of nonlinear growth models in forestry. Silva Fennica 1999;33(4):327-336.
• Firat MZ, Theobald CM, Thompson R. Multivariate analysis of test day milk yields of British Holstein Friesian heifers using Gibbs sampling. Acta Agriculturae Scandinavica, Section A 1997b;47:221-229.
• Firat MZ, Theobald CM, Thompson R. Univariate analysis of test day milk yields of British Holstein Friesian heifers using Gibbs sampling. Acta Agriculturae Scandinavica, Section A 1997a;47:213-220.
• Firat MZ. Bayesian analysis of test day milk yields in an unbalanced mixed model assuming random herd-year-month effects. Turkish Journal of Veterinary and Animal Sciences 2001;25:327-333.
• Forni S, Piles M, Blasco A, Varona L, Oliveira HN, Lobo RB, et al. Comparison of different non-linear functions to describe the Nelore cattle growth. Journal of Animal Science 2009;87:496-506.
• Grimm KJ, Ram N. Nonlinear growth models in Mplus and SAS. Structural Equation Modeling 2009;16(4):676-701.
• Karkach AS. Trajectories and models of individual growth. Demographic Research 2006;15:347-400.
• Kizilkaya K, Balcioglu MS, Yolcu HI, Karabag K, Genc IH. Growth curve analysis using nonlinear mixed model in divergently selected Japanese quails. Archiv Fur Geflugelkunde 2006;70(4):181-186.
• Kizilkaya K, Balcioglu MS, Yolcu HI, Karabag K. The application of exponential method in the analysis of growth curve for Japanese quail. Archiv Fur Geflugelkunde 2005;69(5):193-198.
• Lunn D, Thomas A, Best N, Spiegelhalter DJ. Win BUGS - a Bayesian modeling framework: concepts, structure, and extensibility. Statistics and Computing 2000;10:325-337.
• Marks HL. Growth curve changes associate with long term selection for body weight Japanese quail. Growth 1978;42:129-140.
• Mignon-Grasteau S, Beaumont C, Rochambeau H, Poivey JP, Ricard FH. Genetic parameters of growth curve parameters in male and female chickens. British Poultry Science 1999;40:44-51.
• Narinc D, Aksoy T, Karaman E, Karabag K. Effect of selection applied in the direction of high live weight on growth parameters in Japanese quails. Akdeniz Üniversitesi Ziraat Fakültesi Dergisi 2009;22(2):149-156.
• Narinc D, Karaman E, Firat MZ, Aksoy T. Comparison of non-linear growth models to describe the growth in Japanese quail. Journal of Animal and Veterinary Advances 2010;9(14):1961-1966.
• R Development Core Team. R: a language and environment for statistical computing, reference index version 2.12.1. Vienna: R Foundation for Statistical; 2010.
• SAS Institute. SAS/STAT user's guide. Version 9.2. Cary; 2009.
• Sezer M, Tarhan S. Model parameters of growth curves of three meat-type lines of Japanese quail. Czech Journal of Animal Science 2005;50(1):22-30.
• Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society, Series B 2002;64:583-640.
• Spiegelhalter DJ, Thomas A, Best NG, Lunn D. WinBUGS user manual: Version 1.4. MRC Biostatistics Unit: Cambridge; 2003.
• Theobald CM, Firat MZ, Thompson R. Gibbs sampling, adaptive rejection sampling and robustness to prior specification for a mixed linear model. Genetics Selection Evolution 1997;29:57-72.
• Tzeng RY, Becker WA. Growth patterns body and abdominal fat weights in male broiler chickens. Poultry Science 1981;60:1101-1106.
• Vuori K, Stranden I, Sevon-Aimonen ML, Mantysaari EA. Estimation of non-linear growth models by linearization: a simulation study using a Gompertz function. Genetics Selection Evolution 2006;38:343-358.
• Zhang Z, Hamagami F, Wang L, Nesselroade JR, Grimm KJ. Bayesian analysis of longitudinal data using growth curve models.International Journal of Behavioral Development 2007;31(4):374-383.

Publication Dates

• Publication in this collection
  Jul-Sep 2016

History

• Received
  June 2015
• Accepted
  July 2015