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

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.


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.
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.

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: , ,n 1 2 where y j 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: Richards model: Logistic model: 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).

Bayesian Inference
The following distributions of the three growth model priors ar proposed: Gompertz model: .
Hence, the prior distribution is: where β 0 , β 1 and β 2 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 υ and s 2 ~, s where υ is the degree of belief and s 2 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 υ and s 2 to zero will yield Jeffreys' prior for σ 2 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: (1) Richards 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, σ 2 , was proportional to inverted χ 2 distribution, Suppose we take the Richards model as an example.The full conditionals for the unknown parameters , , , Conditional posterior distributions of the growth curve parameters (β 0, β 1, β 2, β 3 ) 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, σ 2 , 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 nonstandard 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 i θ be the parameters of the i'th model, the DIC is computed as follows (Forni et al., 2009): where, 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.2version (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 (β 0, β 1, β 2, β 3 ), 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 Başar EK, Narinc D

Bayesian Analysis for the Comparison of Nonlinear Regression Model Parameters: an Application to the Growth of Japanese Quail
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.
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).
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 β 0 is higher than its SE in NLS estimation, smaller SEs are obtained for β 1 and β 2 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 β 1 and β 3 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   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, β 3 , allows the flexibility in the inflection point.Moreover, the Gompertz and Logistic growth functions can be derived from the Richards function.When the β 3 parameter in Richards function approaches 0, Gompertz function is approximated, where as if β 3 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.

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

Firat
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, β 0 , 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 β 0 parameter are higher thanAnthony et al. (1991) have reported.Similar results have been reported byAkbas & Oguz (1998),Kizilkaya et al. (2005) andAlkan 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.

Table 2 -
Posterior summary statistics for the three growth functions Sd: Standard deviation; MCError: Monte Carlo Error; dIC: deviance information criterion

Table 1 -
NLS estimation results for the three growth functions AIC: Akaike's information criterion; BIC: Bayesian information criterion