A Bayesian analysis of the effect of selection for growth rate on growth curves in rabbits

Gompertz growth curves were fitted to the data of 137 rabbits from control (C) and selected (S) lines. The animals came from a synthetic rabbit line selected for an increased growth rate. The embryos from generations 3 and 4 were frozen and thawed to be contemporary of rabbits born in generation 10. Group C was the offspring of generations 3 and 4, and group S was the contemporary offspring of generation 10. The animals were weighed individually twice a week during the first four weeks of life, and once a week thereafter, until 20 weeks of age. Subsequently, the males were weighed weekly until 40 weeks of age. The random samples of the posterior distributions of the growth curve parameters were drawn by using Markov Chain Monte Carlo (MCMC) methods. As a consequence of selection, the selected animals were heavier than the C animals throughout the entire growth curve. Adult body weight, estimated as a parameter of the Gompertz curve, was 7% higher in the selected line. The other parameters of the Gompertz curve were scarcely affected by selection. When selected and control growth curves are represented in a metabolic scale, all differences disappear.


INTRODUCTION
Growth curves can describe the entire growth process in terms of a few parameters having a biological interpretation. Selection for growth rate can modify these parameters, but there are some technical difficulties for comparing curves before and after selection. Typically, growth curves are fitted by nonlinear regression or by linear regression if the model can be linearized by transformation (e.g., using a logarithmic scale). The logarithmic scale requires some assumptions: first, the errors are supposed to be multiplicative instead of additive, and, second, it is not possible to find the standard errors of the parameters in the original scale, and approximate standard errors should be used. Moreover, when a Gompertz or a Richards curve is used, a linear form does not exist. When nonlinear regression is used, comparisons between growth curves are not possible because the sampling distribution of the parameters is not known, and approximate methods should also be used. A further difficulty comes from the need to account for possible systematic environmental effects or for genetic relationships between individuals, affecting the structure of the errors. Among the curves proposed, the Gompertz growth curve is widely used to describe the growth of mammals, and it fits better than the other curves for describing the growth of rabbits (Gómez and Blasco [14]). Growth curves have been fitted in rabbits by Baron et al. [2], Fl'ak [8], Rudolph and Sotto [22], Blasco et al. [4] and Blasco and Gómez [5], but only Blasco et al. [4] examined the consequences of selection for growth rate in rabbit growth curves. However, this last study was made without any population control and its results have a limited validity. Some studies draw predictions about the possible correlated response to selection from the heritabilities and correlations (Denise and Brinks [7] in beef cattle; Kachman et al. [15] in mice, Barbato [1] in chickens), but no other studies compare the effect of selection for growth rate on growth curves.
Piles et al. [19] found a positive response to selection in a population of rabbits selected for growth rate. The objective of this research is to examine the effect of selection for an increased growth rate of the rabbit on their growth curve by using a Bayesian procedure derived from the methodology of Varona et al. [26], that overcomes all these difficulties. Other approaches based on linear random regression methods have been suggested (Meyer,[17]), but they are not based on models constructed from the biological meaning of their parameters, as growth curves are. We propose here a nested growth model in which the parameters of the curve are linear functions of environmental and genetic effects. We used a Bayesian inference to assess the correlated response on the growth curve parameters, and the marginal posterior distributions of all unknowns were estimated by Monte Carlo Markov Chain methods. We tested the goodness of fit by using a method that avoids the problems of methods like R-square, strongly dependent on the last part of the curve due to a scale effect. Finally, we expressed the growth curves in Taylor's metabolic scale to better understand how selection for growth rate acts on the live weight growth curve.

Animals
Rabbits come from a synthetic line selected for an increased growth rate. The genetic composition and selection process have been described by Piles et al. [19]. After weaning, rabbits were housed in flat-deck cages, eight rabbits per cage, until they were 9 weeks old, and they were fed ad libitum with a commercial diet (16.0% crude protein, 15.5% fiber, 3.4% fat). Then they were placed in individual cages and the same food was restricted to approx. 140 g per day, since this is the common practice in commercial conditions. At 20 weeks of age they were placed in individual flat-deck reproductive cages, and a commercial diet (17.5% crude protein, 14.5% of fiber and 3.4%) with the same restriction was given.
Embryos from generations 3 and 4 were frozen and thawed to be contemporary of rabbits born in the 10th generation. Offspring from these thawed embryos constituted the control group (C), and were contemporaries to the offspring from parents born in the 10th generation of selection (selected group, S). A total of 137 animals from these groups were individually weighed twice a week the first four weeks and once a week until 20 weeks of age. Males were weighed weekly until 40 weeks of age. The data of the females over 20 weeks of age were not included because they were later pregnant and this modified their growth curves. The numbers of animals measured per group were 27 males and 34 females for group C, and 27 males and 49 females for group S.

Growth model
We describe here a hierarchical model in which each individual i has n i longitudinal data (i.e., the weights from birth to the moment in which the animal died, the individual was eliminated or the experiment stopped). The first stage of the model is the trajectory, and we assumed that the individual growth curve is correctly described using the Gompertz function. The second stage describes how trajectories vary among individuals, and we assumed that growth curve parameters are suitably described by a linear model that includes environmental and genetic effects. A third stage is needed, since a Bayesian probability model requires assigning prior distributions to all unknown quantities.

First stage of the model: the trajectory
We assumed that the weights of each individual follow the Gompertz law: where y ij is the observed weight of the individual i on time j; a i , b i , k i , are the parameters of the Gompertz function for the ith animal, i = 1, 2, . . . , N, and ε ij the residual. Not all individuals have the same amount of records, thus j = 1, 2, . . . , n i . We assumed that the residuals were normally distributed and independent. Other error structures can be proposed; for example, there may be a first-order autoregressive process with heterogeneous variance across the times at which the measurements are taken (Sorensen and Gianola, [24]), and although there is no theoretical difficulty in estimating the parameters in a Bayesian context, this complicates the MCMC process.
We assumed that all animals have the same residual variance at the same time j, but because of a scale effect, the residual variance increases with time until the adult weight is raised, and then remains constant. This can be represented in several ways. After some exploratory analyses fitting the rough data with a Gompertz curve, and examining the s.d. of the residuals, we concluded that the evolution of the standard deviation of the residuals could be represented following a Gompertz law; i.e.: (2)

Second stage of the model: variation among individuals
Each parameter of the curve that describes the trajectory of the growth of each animal is determined by an effect of sex (male or female) and group (C or S), and an environmental component that we assume normally distributed. Calling a, b, k the vectors containing the growth curve parameters where β a , β b and β k are the sex-group effects for the parameters of the growth curve, X is an incidence matrix, and R ⊗ I is the (co)variance matrix of the random environmental effects, where R is the 3×3 (co)variance matrix between the residuals e a , e b , e k , and I is a N × N identity matrix. This means that, for each individual i: This assumption is based on the biological meaning of the parameters: if a describes the adult weight and k is related to the slope of the curve (the growth rate), it is reasonable to suppose that they will be correlated within individuals, but not between individuals, given that the genetic relationships between individuals are not considered at this stage. We simplify now the notation, naming β = [β a , β b , β k ] the vector with all sex-group effects, and p = [a , b , k ] the vector with the Gompertz curve parameters for each animal. We name p ε = [a ε , b ε , k ε ] the vector with the Gompertz curve parameters for the s.d. of the residuals, thus: and we will call this model, Model 1.

Third stage of the model: uncertainty about the second stage parameters
We consider that the sex and group effects have a normal prior distribution: where m and V are the subjective mean and variance for the prior beliefs about the systematic effects. We propose, according to Sorensen et al. [23], an inverted Wishart distribution for prior distributions of R: where (S R , n R ) are the hyperparameters of the inverted Wishart function. These hyperparameters, modify the shape of the function changing the amount of information of the prior density (see Blasco [3] for a detailed discussion about the prior information). The prior distributions for the parameters a 0 , b 0 , k 0 of the residual standard deviation are assumed to be flat with limits that guarantee the property of the distribution. We always used proper prior distributions in order to guarantee all the posterior distributions to be proper. Model 1 ignores that data are correlated because it does not take into account the genetic relationships between individuals. This produces an underestimation of the standard deviation of the posterior densities. A model including all the genetic effects of all animals from the first generation of selection has been proposed by Varona et al. [27] for dairy cattle. We cannot use this model here because we only have data for growth curves from the last generation of selection and from generations 3-4. In order to assess the effect of the relationships between animals, we fit a model in which the growth curve parameters of each individual were also determined by a genetic effect. We will call this model, Model 2: (p|β, u, R) ∼ N(Xβ + Zu, R ⊗ I); (5) where the genetic effects have a normal prior distribution: (u|G) ∼ N(0, G ⊗ A); (6) where G is the genetic (co)variance matrix between the Gompertz growth curve parameters and A is the numerator relationship matrix including only the relationships of the individuals of groups C and S. The prior distribution of G is also an inverted Wishart distribution:

Bayesian inference
The joint posterior distribution is (Appendix A) Prior distributions represent the state of knowledge before the results of the experiment become available. For the group effects β we have used vague priors, taking m and V from a previous experiment of Blasco and Gómez [5], who estimated the growth curve of this line in the base generation. Since there is no information on the residual (co)variances for the growth curve in rabbits, two different priors were used to express a vague knowledge about the (co)variance matrix R. We can then compare the two possible states of opinion, and study how the use of the different prior distributions affects the conclusions from the experiment. We first used flat priors (with limits that guarantee the property of the distribution) for two reasons: to show an indifference about their value and to use them as reference priors, since they are usual in Bayesian analyses. Since prior opinions are difficult to draw in the multivariate case, we chose the second prior by substituting a (co)variance matrix of the components in the hyper parameters S R and S G and using n R = n G = 3, as proposed by Gelman et al. [11] in order to have a vague prior information. These last priors are based on the idea that S is a scale-parameter of the inverted Wishart function, thus using for S R and S G prior covariance matrixes with a low value for n, would be a way of expressing prior uncertainty. We proposed S R and S G from phenotypic covariances obtained from the data of Blasco and Gómez [5]. Table I shows the hyper parameters of both prior distributions.
The conditional distributions needed to run the Gibbs sampler are derived in Appendix B. Conditional posterior distributions for a i and β i and u are Normal distributions, conditional posterior distributions for the (co)variance components (R and G) are Inverted Wishart distributions and conditional posterior distributions of b i , k i and p ε are non standard statistical distributions. There are algorithms for the exact random sampling of Normal and Inverted Wishart distributions, but when the distribution is not a standard one, an  i and k (t) i [11]. After several trials, the proposal distributions used were: where U is the uniform distribution. The choice of the limits for that distribution determines the acceptance rate. If the width of such an interval is too small, the proposed values will be closed to the current ones, the rejection rate will be low but the process will move slowly throughout the parameter space. On the contrary, if it is too large, the proposed values are far away from the current ones and this results in a high rejection rate. The scale of the proposal distribution was determined in a preliminary analysis. The above choice led to acceptance rates ranging between 17 and 45%. For each analysis three chains with different starting values were run. After several trials, the length of each chain was set to 300 000. The burn-in period was 150 000 iterations, higher than the minimum burn-in required according to the method of Raftery and Lewis [20], and the sampling interval was 10, so that a total of 15 000 samples were kept from each chain. Convergence was tested for each chain separately using the criterion of Geweke [12]. Convergence was also assessed by the test of Gelman and Rubin [10]. For each variable, a scale parameter ( √ R), which involves the variance between and within the chain was computed. This parameter can be interpreted as the factor by which the scale of the marginal posterior distribution of each variable would be reduced if the chain were run to infinity, and it should be close to 1 to convey convergence. All these samples were used to estimate the features of posterior distributions. Autocorrelation between samples and Monte-Carlo error of features of marginal distributions [13] were also calculated.

Outliers
Preliminary analyses were conducted to detect the presence of outliers or atypical growth patterns. An observed weight was declared to be an outlier if the standardized absolute value of the residual posterior mean was larger than three standard deviations from the standard normal distribution [6]. An atypical growth pattern was declared when the Mahalanobis distance between the individual growth curve parameters and the average of its group was high. Since we have three parameters, the square of this distance 3 . We checked how many individuals had a value of D 2 laying in the area of P < 0.01.

Goodness of fit
The goodness of fit was checked by the square of the correlation coefficient between the predicted and observed values r[E(Y r |y −r ), y r ]. This global criterion, like the coefficient of determination of the fit, has the disadvantage of depending more on the last part of the curve than on the first part due to a scale problem because the absolute value of the errors are higher at the adult weight than at the beginning of growth. Moreover, nonlinear models require to examine the whole growth trajectory, since a growth curve can fit well in some parts but not in others. Due to this, we used cross-validation predictive densities to asses the goodness of fit of the model. The observed values y r were compared with their prediction Y r obtained using all the other data y −r . We used one of the checking functions proposed by Gelfand et al. [9]: We obtained E(g|y −r ) for each observed value r. This expectation shows the probability of a predicted value of being higher or lower than the observed one. If the model fits the data properly, E(g|y −r ) should be close to 0.5, thus a global criterion for goodness of fit is to calculate the average of these expectations for all individuals in each point t j of the growth curve. A graph with these averages shows whether the fit is good along the curve or whether there are parts of the curve that fit better than others. This technique has the advantage of being free of the scale effect. The expectation of the checking function E(g|y −r ) is computed using the MCMC methods. These methods are very computing demanding, thus we applied an importance sampling procedure as suggested by Rekaya et al. [21].  Figure 1 shows the weekly averages of live weight of males (M) and females (F) of the control (C) and selected (S) groups. Nine observations were declared to be outliers and were removed. No animals presented an atypical growth pattern. Group S showed a higher live weight for both males and females along the whole growth curve. No sexual dimorphism was observed at slaughter time (9 weeks of age), but this dimorphism appeared at later ages, females being heavier than males, as in rabbits and not in other domestic species.

RESULTS
The analyses of growth curves made with two different priors gave very similar results, showing that the information of these analyses come essentially from the data and not from the priors used. Since the results from both priors were almost the same, only the results obtained using the flat prior will be commented. Tables II and III show the means and standard deviations of the posterior densities of the curve parameters for the flat prior, as well as the Monte Carlo standard errors and convergence tests of the Gibbs sampler for the growth curves. The autocorrelation was generally low, in the model without genetic effects, but it was higher in the model with genetic effects, leading to higher estimates of Monte Carlo standard errors. All chains gave very similar results, the difference between chains being of the same size of the Monte Carlo Standard error, thus they were blended to give the estimates of the means and s.d. The convergence was good, the z-score of the Geweke test in the model without genetic effects was generally low, and never higher than 1.96, and the scale parameter of the Gelman and Rubin test was always close to 1. The model  with genetic effects showed one case in which the z-score was higher than 1.96, but the results of the Gelman and Rubin test were good. The burn-in period used was much higher than the minimum recommended by the procedure of Raftery and Lewis. Thus no pathologies were detected in the sampling process. The square of the correlation coefficient between the predicted and observed values was 0.99. Figure 2 shows the averages of the expectations of the Gelfand checking function for each point of the growth curve. Although all of them lay  near 0.5 (none of the predicted values showed a high probability of being lower or higher than the observed value), it can be seen that at the beginning there is a trend of obtaining predictions higher than the observations, and at the end of the growth curve there is a certain trend of obtaining predictions that are lower than the observations. Adult weight is raised very slowly and the asymptotes of the growth curves tend to slightly underestimate the final adult weight.  Figure 3 shows the posterior distributions for the differences between the curve parameters of lines S and C. All the distributions were approximately normal, similar to the ones shown in this figure. The parameters for males were estimated with a higher accuracy, since they have more individual data. Tables II and III also show the differences between growth curve parameters of groups S and C for males and females. Although the estimated means generally agree, standard deviations were, as before, higher for the model including the genetic effects. Group S showed a higher parameter a but the differences found for parameters k and b were very small. The effect of selection seems to be a scale effect, increasing all weights along the growth curve but not changing the shape of the curve, as predicted by Taylor [25]. Figure 4 shows the fitted growth curves for males in their original scale and in the metabolic scale of Taylor [25]. All the selection effects disappear when the curves are represented in Taylor's metabolic scale.

DISCUSSION
The Bayesian inference makes use of prior information. In the multivariate case, it is almost impossible to define with some accuracy this previous knowledge, because of the difficulty of expressing prior opinions in a multidimensional space (see Blasco, [3] for a detailed discussion about prior information in a Bayesian context). We have compared here two prior opinions for the dispersion parameters, one expressing indifference, and the other based on prior univariate information derived from the experiments from another data base of the same breed. Although both were constructed from different bases, they    give virtually the same results, which means that the experiment had enough data and the prior information was rather irrelevant to the results obtained. Monte Carlo Markov chain techniques have solved the difficult operative problems that have prevented the application of Bayesian techniques in genetics for many years. One of the advantages of the use of these chains is that new samples of marginal posterior distributions of new variables can be derived from the chains obtained in the Gibbs sampling process. Thus, to compare parameters from two different growth curves, the marginal posterior density of the difference between parameters can be obtained just by calculating the difference between chains in each sample.
Model 1 has less information than model 2, but it does not depend on the genetic estimation of the variance components, which are necessarily inaccurate given the reduced amount of data common in this type of experiment. Since genetic effects are not included in model 1, this tends to underestimate the s.d. of the marginal posterior distribution. Model 2 has the advantage of considering the relationships between individuals, however this model depends on the estimates of the genetic parameters. These estimates cannot be very accurate given the limited amount of data, and there is very little in the literature that can help (few experiments, also with few data). Including the whole relationship matrix and having data of individuals in all generations, the correlated response to selection for growth rate on the curve parameters could have been estimated without the need of a control population. However, we only had data in the last generation of selection, not in the other ones, and the lack of data in all the other generations makes this task particularly difficult. Due to this, we used a relationship matrix that only included the individuals of groups C and S and their parents. This takes into account the main relationships, which permits to calculate more accurately the s.d. of the marginal posterior distributions, but avoids to estimate all parameters of all individuals from the first generation of selection. The results of model 1 and model 2 are remarkably similar, showing that including or not the genetic effect has a small effect in the estimation of the means of the marginal posterior distributions.
No other authors have studied hitherto the effect of selection for growth rate on the growth curve in rabbits. Blasco et al. [4] fitted rabbit growth curves in an unselected population, and in the same population ten generations later, but the absence of a control population makes their results merely indicative. Even in other species, the studies of the effect of selection on the growth curve has been approached only indirectly, based on the estimation of the genetic parameters of the growth curve more than in the direct comparisons of the effect of selection in these parameters [15,18]. We have exemplified here a way of comparing growth curves between populations by using a procedure to assess the differences between the growth curve parameters.
An undesirable consequence of selection for growth rate is the increment of the adult weight because it augments the costs of maintaining a parent population. Adult weight (a-parameter of the Gompertz curve) increased with selection, whereas the parameters related to the slope of the curve did not practically change (Fig. 4). This result is more clear when comparing males of group S and C than when comparing females, probably due to the limited amount of records in the females side. When curves are represented in the Taylor's metabolic scale all the effect of selection disappears (Fig. 4). This was predicted by Taylor [25], who stressed that all weights are correlated and a selection for growth rate can easily lead to an increase in adult weight with no changes in the curve slope. Changes in the curve slope have been produced by McCarthy and Baker [16] in mice, but the selection process was very inefficient, since the remaining genetic variation after restricting adult weight is very small due to the genetic correlations between all weights along the growth curve. In this circumstance, it can be predicted that male lines will become giant lines and the management of reproduction will be more difficult, unless artificial insemination is used.

APPENDIX A
To construct the posterior distribution we need, according to Bayes Theorem, the likelihood and the prior distribution. Since we have considered the residuals ε ij independent, the likelihood is a product of normal functions: f (y|p, β, u, G, R, p The notation of the likelihood can be simplified, since given the parameters p implies that β, u, G and R are also given. For the prior distributions, we first consider that prior for p ε is independent from the other priors. Then we consider that G and R are independent, and β and u are also independent. We also take into account that given u implies that G is given and the notation can be simplified. The posterior distribution is: f (p, β, u, G, R, p ε |y) = f (y|p, p ε )f (p, β, u, G, R)f (p ε )/f (y) = f (y|p, p ε )f (p|β, u, R)f (β)f (u|G)f (G)f (R)f (p ε )/f (y).
We can put this function in an explicit form, since f (G) and f (R) are Inverted Wishart distributions, f (p ε ) is a uniform distribution, and the other functions are normal distributions. The likelihood is, from (1) where σ j has the expression (2) function of the components of p ε . Calling the likelihood of the data of individual i is: where Σ i is the (co)variance matrix corresponding to individual i. The prior distributions are, from (3), (4), (5), (6) and (7): where g −i is the vector g = [β , u ] with the ith element excluded, RHS i is the corresponding element of the Right Hand Side and c ij the corresponding coefficient of the mixed model equations constructed as if the observed traits were the Gompertz curve parameters: To access this journal online: www.edpsciences.org