Estimates of direct and maternal covariance functions for growth of Australian beef calves from birth to weaning

Records for birth and subsequent, monthly weights until weaning on beef calves of two breeds in a selection experiment were analysed fitting random regression models. Independent variables were orthogonal (Legendre) polynomials of age at weighing in days. Orders of polynomial fit up to 6 were considered. Analyses were carried out fitting sets of random regression coefficients due to animals' direct and maternal, additive genetic and permanent environmental effects, with changes in variances due to temporary environmental effects modelled through a variance function, estimating up to 67 parameters. Results identified similar patterns of variation for both breeds, with maternal effects considerably more important in purebred Polled Herefords than a four-breed synthetic, the so-called Wokalups. Conversely, repeatabilities were higher for the latter. For both breeds, heritabilities decreased after birth, being lowest when maternal effects were most important around 100 days of age. Estimates at birth and weaning were consistent with previous, univariate results.


INTRODUCTION
Random regression (RR) models and the associated covariance functions (CF) have been advocated for the analysis of "traits" measured repeatedly per individual. They facilitate modelling changes in the trait under consideration over time and in its (co)variance structure. Applications so far have concentrated on the analysis of test day records of dairy cows [4,8,[12][13][14]30,31,34,35,37,39,40]. Other work considered growth and feed intake in pigs [3,36], feed intake and live weights in dairy cattle [41] and weights of mature beef cows [24,25].
Previous applications fitted at most two sets of RR coefficients per animal, namely due to animals' direct genetic and permanent environmental effects. Early growth of mammals, however, is subject to substantial maternal effects in addition, both genetic and environmental. Estimation of maternal effects has been found inherently problematic, even for "standard", univariate analyses of traits such as birth or weaning weight. Extension of RR models to include maternal effects is conceptually straightforward. However, inclusion of additional sets of RR coefficients to accommodate maternal effects increases the complexity of corresponding analyses considerably, especially if we want to distinguish between genetic and permanent environmental maternal effects. This paper presents a RR analysis of weights of beef calves from birth to just after weaning, attempting to separate direct and maternal covariance functions.

Data
Records originated from a selection experiment in beef cattle carried out at the Wokalup research station in Western Australia. This experiment comprised two herds of about 300 cows each. The first were purebred Polled Hereford (PH), and the other a four-breed synthetic formed by mating Charolais × Brahman bulls to Friesian × Angus or Friesian × Hereford cows, the so-called "Wokalups" (WOK), see Meyer et al. [26] for details. Management of the experiment comprised short mating periods (natural service) of 7 to 8 weeks, resulting in the bulk of calves being born during April and May each year. Calves were weighed at birth and subsequently, from July to late December or early January, at monthly intervals. This resulted in up to nine weight records per calf. Calves were weaned between mid-November and early January each year, i.e. the last weight represented a post-weaning weight in most years. Variation in weaning dates was necessitated by seasonal conditions -the data included several years of drought. Consequently, annual means in age at weaning ranged from 182 to 246 days, with overall means of 211 and 214 days for PH and WOK, respectively.
Data selected consisted of birth and subsequent, monthly weights for calves born between 1975 and 1990. This yielded 21 272 and 22 230 weights for PH and WOK, respectively. Basic edits eliminated records with implausible dates or weights. In addition, changes in weights between individual weighings were scrutinised. The number of animals in the data was sufficiently small to allow questionable sequences of records to be inspected individually. Apparently aberrant records, in particular records clearly out of sequence, identified on basis of both average daily gain as a proportion of the mean weight and absolute change in weight, were disregarded. In doing so, allowance was made for large variation in gain between birth and the first post-natal weight and an increased chance for a decrease in weights between November and December recordings due to weaning stress. Limits were established based on means and distributions of average daily gains for each monthly weighing (across years). Changes in average daily gains as proportion of the mean as large as −1.3% to 3.9% between birth and first weighing and −0.8% to 1.7% between weights preand postweaning were allowed, while corresponding changes between other monthly weighings were restricted to −0.3%-0% and 1.2% to 2.0%. Figure 1 shows the distribution of records over ages at weighing. Ages ranged up to 297 days. However, there were few records after 280 days of age. These were thus eliminated to avoid problems due to small numbers of records per subclass. This left 21 053 and 21 807 records for PH and WOK, respectively, on 3 417 (PH) and 3 768 (WOK) animals. Birth weights were available for almost all calves (3 406 for PH and 3 727 for WOK). Further characteristics of the data structure are given in Table I.

Preliminary analyses
"Standard" univariate, animal model analyses were carried out to assess changes in variance components due to direct and maternal effects with age. These considered single records per animal only. Records were selected for target ages at fortnightly intervals with a separate class for birth weights. This yielded 19 partially overlapping data sets with ages of 0, 2-35, 21-49, 35-63, . . . 245-259 and 245-280 days.
Two analyses were carried out for each data set. The first fitted a model with permanent environmental maternal effects in addition to animals' additive genetic effects only (Model U1). The second model included both genetic and environmental maternal effects (Model U2). Corresponding variance components were estimated by REML. Differences in likelihoods for each pair of analyses were used to assess the importance of maternal genetic effects, and the ability to separate environmental and genetic maternal components. Fixed effects for univariate analyses were as for RR analyses (see below), with a linear regression on age at weighing fitted within sex in addition.

Fixed effects
Mean age trends were taken into account by a fixed, cubic regression on orthogonal polynomials of age (in days). Preliminary investigations had shown higher orders of fit to yield virtually no reduction in residual sums of squares, presumably due to a close association between age at recording and contemporary group subclasses. The same order of fit for the fixed regression on age was considered throughout, making REML likelihoods for different orders of fit of RR directly comparable.
Other fixed effects fitted were similar to those fitted in earlier, non-RR analyses of birth, weaning and later weights for these data [26]. These included contemporary groups (CG), defined as paddock-sex of calf-yearweighing number (1 to 9) subclasses and a birth type (single vs. twin) effect. Dam age was modelled as a yearly age class effect. As in previous analyses, Wokalups were treated as a breed and no specific crossbreeding effects were fitted.

Random effects
All RR models fitted Legendre polynomials (e.g. [1]) of age at recording (in days) as independent variables. Orders of polynomial fit up to k = 6 were considered. Polynomials included a scalar term, i.e. involved powers of age up to five. RR analyses fitted a set of k regression coefficients for each random effect considered, up to four in total. For simplicity, k was initially chosen to be the same for all factors.
The first set of analyses fitted three sets of RR coefficients, namely due to animals' direct genetic effects (A), due to animals' permanent environmental effects (R) and due to dams' permanent environmental effects (C) (Model G1). The second set fitted a fourth set of RR coefficients in addition, namely due to maternal genetic effects (M) (Model G2). Both models incorporated all pedigree information available, assuming direct and maternal genetic effects (A and M) were distributed proportionally to the numerator relationship matrix between animals.
When comparing orders of fit of polynomial equations, it is recommended to consider the next higher order of fit involving the same type of exponent (odd versus even) [10], p. 182-183. For instance, if a cubic polynomial fitted the data, the linear coefficient was likely to explain a significant amount of variation while the quadratic term might contribute nothing. Hence k was initially increased in steps of 2, starting at a cubic regression (k = 4), as preliminary analyses showed linear (k = 2) and quadratic (k = 3) regressions to be clearly inadequate to model the variation in the data [23]. Orders of fit greater than k = 6 were not examined as preliminary work had also indicated that this would be unnecessarily high.
Further analyses considering different orders of fit for the four random effects where carried out subsequently, with the choice of models determined by results of the earlier analyses. The aim in doing so was to determine the minimum order of fit required for each random factor, and thus determine the most parsimonious model describing the data.

Additional analyses
Results from the random regression analyses raised questions about the influence of postweaning weights, as well as records at the highest ages and birth weights on the minimum orders of fit required. Hence, additional analyses were carried out for PH, considering subsets of the data. First, as weaning dates were known, all weights taken post weaning were eliminated. In addition, weights at late ages (> 250 days) were discarded and animals with birth weight records only which were deemed to contribute little information were omitted. This left 19 399 records on 2 865 animals, with 3 260 animals in the analysis and 2 014 fixed effects levels (rank 2 011) in the mixed model equations. Secondly, birth weight records and records up to 14 days of age were disregarded in addition, reducing the data to 16 438 records on 2 863 animals, with 3 260 animals and 1 727 fixed effects (rank 1723) in the analysis.

Model of analysis
More formally, this gave a model of analysis for y ij , the j-the record on animal i of with a * ij denoting the age at recording for y ij , standardised to the interval [−1 : 1] for which Legendre polynomials are defined, and φ m (a * ij ) the corresponding m-th Legendre polynomial. F ij represented the fixed effects pertaining to y ij , and b m the coefficients of the fixed, cubic regression modelling mean age trends. β Qm was the m-th random regression coefficient for the Q-th random effect, with Q standing in turn for A, R and C for model G1, and for A, M, R and C for model G2, and k Q the corresponding order of polynomial fit. Finally, ε ij denoted the residual error.

Variance and covariance functions
Parameters estimated in RR analyses were the matrices of covariances between RR coefficients: Covariances between RR coefficients pertaining to different random factors were assumed to be zero throughout.
Elements of K Q are the coefficients of the covariance function defining covariances between any two ages in the data for the Q-th random effect (e.g. [28]). Hence, estimated covariances between records for animal i at ages a ij and a ij are with K Q mn denoting the mn-th element of K Q , and Ω Q = ω Q mn defining the Q-th covariance function [16,22]. In addition, (1) allowed for temporary environmental effects or "measurement errors", ε ij . These were considered to be independently distributed throughout with variances σ 2 . Commonly it is assumed that these variances are homogeneous, consistent with the concept of a true measurement error affecting all observations equally. Preliminary analyses, however, had clearly shown that this assumption as inadequate [23]. Hence σ 2 was considered to change with age, with changes described by a polynomial variance function (VF) with σ 2 j the variance at the j-th age, σ 2 0 denoting the measurement error variance at the mean age (0 on the standardised scale) and b r the coefficients of the VF. The VF has v + 1 parameters, comprising the v coefficients b r and σ 2 0 . An alternative VF entails regression of the logarithm of the variance on age (e.g. [7,11,33]). This of particular interest as it allows an exponential increase in variance with time to be modelled parsimoniously by a simple linear regression. Moreover, in contrast to (4), it does not require any constraints on the coefficients of the polynomial regression to be imposed.

Estimation
Estimates were obtained by restricted maximum likelihood (REML) using program "DXMRR" [21], employing a combination of average information (AI-REML) and derivative-free algorithms to locate the maximum of the likelihood. REML estimation of covariances between RR coefficients is analogous to multivariate estimation in "standard" (i.e. non-RRM) analyses [22].
AI-REML algorithms for the latter have been described by Madsen et al. [18] and Meyer [20]. Additional calculations required to estimate the parameters of a VF for measurement error variances with an AI-REML algorithm are outlined in the Appendix.
Search for the maximum of the likelihood was invariably slow. With highly correlated parameters, several restarts were required for each analysis. The AI-REML algorithm tended to perform less well than in equally dimensioned, non-RR multivariate analyses. If estimates of covariance matrices had eigenvalues less than 0.001 these were set to an operational zero (10 −7 ) and estimation was continued fixing these values, effectively forcing estimated matrices to have correspondingly reduced rank (r) (see [22]). Generally, this resulted in improved convergence of the iterative estimation procedure.

Model selection
Fit of different models was compared by examining estimated variances (σ 2 A : direct, additive genetic variance, σ 2 M : maternal, additive genetic variance, σ 2 R : direct, permanent environmental variance, and σ 2 C : maternal, permanent environmental variance) for ages in the data and comparing maximum likelihoods and information criteria for each analysis. To account for nonstandard conditions at the boundary of the parameter space [38] in carrying out likelihood ratio tests (LRT), differences in log L were contrasted to χ 2 values corresponding to twice the error probability of α = 5%.
Information criteria comprised the REML forms of Akaike's Information Criterion (AIC) and Schwarz' Bayesian Information Criterion (BIC). Let p denote the number of parameters estimated, N the sample size, r(X) the rank of the coefficient matrix of fixed effects in the model of analysis, and log L be the REML maximum log likelihood. The information criteria are then given as [43] AIC = −2log L + 2p (6) and

RESULTS
Numbers of records and means for individual ages (weekly intervals) are shown in Figure 1. Almost all animals had birth weight records (not shown). Growth for both breeds was approximately linear. While there was little difference in size at birth, WOK calves grew faster throughout than PH with means (SD) of 157.5 (88.1) and 138.6 (79.1) kg, respectively. Corresponding SD are shown in Figure 2. Values for both breeds were again very similar and increased steadily with age, both on the observed scale and for data adjusted for least-squares estimates of fixed effects. The latter represents the pattern of variation to be modelled by the estimated covariance functions. Corresponding coefficients of variation (CV), decreased with age, i.e. variances increased less than might be anticipated due to scale effects. Except for the highest ages, CVs were thus consistently higher for PH than WOK. Due to computational requirements, only limited comparisons between different orders of fit could be carried out. As suggested by results from preliminary analyses [23], only RR due coefficients due to permanent environmental maternal effects was fitted initially (Model G1). Results from "standard" (i.e. non-RR) analyses of birth and weaning weights in beef cattle, comparing different models of analyses, often showed that, when fitting only one of the maternal effects, this is likely to "pick up" most of the total maternal variation, i.e. due to both genetic and permanent environmental effects (e.g. [19]). It was assumed that the same pattern in partitioning of variation would apply for RRM analyses.
For both breeds, a model fitting Legendre polynomials to order k = 6 for all three random effects, denoted by 6066 in the following (the four digits corresponding to the orders of fit for A, M, R and C, respectively), and a cubic VF for measurement error variances (v = 3) -involving a total of 67 parameters to be estimated -was considered more than adequate on the basis of earlier results [23].

K. Meyer
Estimated covariances among regression coefficients from these analyses (k = 6066, v = 3) showed little variation in the quartic and quintic regression coefficients for direct genetic and even less for maternal, environmental effects. Analyses were thus repeated reducing the order of fit for these effects to a cubic regression, while still fitting a quintic regression for direct, permanent environmental effects (i.e. k = 4064). BIC (see Tab. II) for both breeds were smaller for this model, i.e. suggested that 45 (rather than 67 parameters for k = 6066) sufficed to describe the pattern of variation in the data.
A number of additional analyses involving different combinations of cubic and quintic polynomial regressions for the different random effects (k = 4044, 4064 and 6064) were carried out subsequently. In addition, analyses fitting model G2 considered orders of fit k = 4444, 4464 and 6464. In doing so, choices of k were guided by results from analyses carried out so far, and not all models were fitted for both data sets.
Values for log L and corresponding information criteria for the different analyses are given in Table II. As for univariate analyses, fitting maternal genetic effects for WOK did not increase log L significantly (k = 4464 vs. k = 4064). Values for log L for WOK were highest for the model involving most parameters (k = 6066) though not significantly higher than for k = 6064. The information criteria too suggested that a cubic regression for maternal environmental effects sufficed (k = 6064).
In contrast, fitting model G2 instead of G1 for PH yielded a marked increase in log L, consistent with results from univariate analyses. LRT and both information criteria suggested that of the models examined so far, a quintic regression for both direct effects and a cubic regression for both maternal effects (k = 6464) fitted best for PH. Results from this analysis found little variation in the cubic regression coefficient for maternal, permanent environmental effects. Reducing the order of fit accordingly, i.e. to k = 6463 did not reduce the likelihood significantly.
Further analyses eliminated the highest order RR coefficient for random factors with comparatively small variances (around 1). Whilst this tended to cause a significant decrease in log L, corresponding BIC values decreased which suggested that the reduced model was adequate and provided a more parsimonious representation of the covariance structure in the data. As shown in Table II, log L and the AIC favoured orders of fit of k = 6463 for PH and k = 6064 for WOK, with 62 and 56 parameters, respectively.
It follows from (6) and (7) that BIC imposes a much more stringent penalty for the number of parameters fitted than AIC. For our data, the factor log N − r(X) was close to 10 (9.85 for PH and 9.88 for WOK), i.e. almost five-fold that for AIC. This let models with k = 5163 and 47 parameters and k = 5062 with 43 parameters to be selected as "best" for PH and WOK, respectively.   Reducing the order of fit for direct additive genetic effects to a cubic regression (k = 4263 vs. k = 5263 for PH, k = 4062 vs. k = 5062 for WOK) or eliminating the quintic regression coefficient for direct, permanent environmental effects (k = 5253 vs. k = 5263 for PH, k = 5052 vs. k = 5062 for WOK), however, clearly resulted in a less appropriate model. Similarly, decreasing the order of fit for maternal, permanent environmental effects by one (k = 5262 vs. k = 5263 for PH, k = 5061 vs. k = 5062 for WOK), did not provide sufficient scope to model changes in variation with time any longer, resulting in substantially higher information criteria. Reducing the order of fit for direct effects by two often resulted in a somewhat higher BIC than a reduction of one (k = 5253 and k = 5243 vs. k = 5263 for PH, k = 4263 and k = 3263 vs. k = 5263 for PH, k = 4062 and k = 3062 vs. k = 5062 for WOK, but not k = 5052 and k = 5042 vs. k = 5062 for WOK), suggesting that the cubic regression coefficients for direct genetic effects and the quartic coefficient for direct, permanent environmental effects in PH were of lesser importance.
Whilst BIC was smallest for k = 5163 for PH, this model implied a constant maternal genetic variance (of 4.1 kg 2 ) for all ages. With phenotypic variances increasing with ages, this would yield a maternal genetic heritability with maximum value at birth and decreasing steadily with age. Clearly, this "choice" was due to the stringent penalty for the number of parameters imposed, yielding an only marginally higher BIC if maternal genetic effects were omitted from the model altogether (k = 5063). It emphasized that the data did not support an accurate partitioning of maternal effects into their genetic and environmental components, in spite of the fact that it contained a sizeable number of records for calves born after embryo transfer [26], which were expected to reduce sampling correlations between the two maternal effects. Such patterns of variation as suggested by models k = 5163 or k = 5063, however, clearly disagreed with univariate analyses and or expectations, based on the multitude of analyses reporting estimates of maternal genetic variances for birth and weaning weights available in the literature. Hence model k = 5263 with slightly higher BIC and 49 parameters was considered most appropriate for PH, and treated as "best" in the following.
Cubic variance functions were selected to model changes in measurement error variances with time, again based on preliminary analyses [23]. Initial comparisons between a model fitting a VF for log(σ 2 ) rather than σ 2 showed the latter to be advantageous (higher log L). Hence further analyses were carried out fitting a cubic VF for σ 2 .

Covariance functions
Estimates of covariance matrices between RR coefficients (K Q for Q = A, M, R, C) and corresponding correlations are summarised in Table III for k = 5263 for PH and k = 5062 for WOK. The resulting covariance functions (Ω Q , see (3)) are given in Table IV. In all cases, intercepts of the polynomial regressions were most variable, and there were strong positive correlations between the intercept and the linear coefficient. Conversely, correlations between the intercept and the quadratic regression coefficient were negative except for K A for WOK, ranging from moderate to high.
Strong correlations between RR coefficients caused one eigenvalue of estimated covariance matrices for direct, permanent environmental effects (K R ) to be essentially zero. Fitting both maternal effects for PH resulted in an estimate of K M with one negligible eigenvalue. The first eigenvalue for all CF dominated throughout, indicating that a large proportion of the total variation is explained by the first eigenfunction of each CF.

Genetic parameters
Estimates of variance components for the ages in the data are shown in Figure 3, for both the best model as determined by BIC and by AIC and log L. Corresponding estimates of genetic parameters are given in Figure 4.
Estimates of variance components followed a similar pattern for both breeds, but differed in magnitude. Direct components were consistently higher for WOK than for PH. There was little difference in values from the best model (k = 5263 for PH and k = 5062 for WOK) and the model selected using a LRT or AIC (k = 6463 and k = 6064, respectively). Differences in estimates between the two models were largest for σ 2 M in PH, where the order of fit was reduced from a cubic to a linear regression.
Estimates of σ 2 A , σ 2 M and σ 2 C showed reasonable agreement with corresponding univariate estimates, except for σ 2 A in PH. For both breeds, estimates of σ 2 P from RR model analyses were consistently higher than their counterparts from univariate analyses, increasingly so with increasing age. This discrepancy was slightly larger for PH than WOK, presumably due to higher estimates of σ 2 A . Heritability. Estimates of σ 2 A increased less than σ 2 P for the first 120 days, resulting in a decrease in estimates of the direct heritability (h 2 ) after birth, with a minimum for both breeds at about 100 days of age. For WOK, estimates of h 2 showed good agreement with their univariate counterparts. For PH, however, estimates from RR analyses were consistently higher except at birth, due to Orders of fit for direct additive genetic, maternal additive genetic, direct permanent environmental and maternal permanent environmental polynomials, respectively. Orders of fit for direct additive genetic, maternal additive genetic, direct permanent environmental and maternal permanent environmental polynomials, respectively.  higher estimates of σ 2 A . In part this might be explicable by increased sampling variation in the more detailed models (G2) fitted. When comparing estimates of h 2 for PH from models G1 (not shown), results from the different types of analyses agreed almost as well as for WOK.

K. Meyer
Repeatability. Direct, permanent environmental effects ( p 2 ) explained only about 12% of variance in birth weights, but increased steeply in the first month, reaching 32% (PH) and 45% (WOK) at 31 days of age. For WOK, estimates of p 2 increased further, reaching 53% around 3 months of age, while estimates for PH increased slowly, reaching a value of 38% around 8 months. The increase and subsequent decline in estimates of p 2 for WOK between about 2 and 6 months of age was absent in PH, suggesting that sampling variation in the partitioning of the between animal variation into its genetic and environmental components may have contributed to overestimate of σ 2 A and h 2 . Estimates of the repeatability (t = h 2 + p 2 , not shown) exhibited markedly less changes with age than h 2 and p 2 . This is again indicative of a negative sampling correlation between direct genetic and permanent environmental effects. For both breeds, t peaked at about one months of age (65% for PH, 82% for WOK), declined to a minimum of 57% at about 4 months for PH and 77% at 5 months for WOK, and then increased again to 62% (PH) and 80% (WOK) at 8 months of age.
Maternal effects. Maternal effects were considerably more important in PH than in WOK. Estimates of σ 2 C for PH were consistently larger than of σ 2 M . Values for c 2 = σ 2 C /σ 2 P increased steadily from birth to about 3 (PH) to 4 (WOK) months of age, being highest when direct genetic effects were least important, and decreasing slowly from there. Small increases in estimates after day 250 were accompanied by corresponding declines in estimates of σ 2 P , presumably reflecting reduced variation for these ages with few records (cf. Fig. 1). Estimates of the maternal genetic heritability (m 2 ) for PH varied little over the range of ages considered and showed reasonable agreement with estimates from univariate analyses. As for c 2 , an increase in estimates after 250 days of age was considered to be an artifact of a decrease in σ 2 P due to small numbers of records.

Correlations
Estimates of direct genetic (r A ), direct permanent environmental (r R ), maternal genetic (r M ), maternal permanent environmental (r C ), and phenotypic (r P ) correlations among the ages in the data for PH (k = 5263) are shown in Figure 5. Whilst variance components for individual effects differed markedly for the two breeds, patterns of correlations were very similar. Taking correlations as depicted in Figure 5 (2-day intervals from 0 to 10 days, 5 day intervals from 10 to 280 days) and considering the lower triangle of the correlation matrix only  Genetic Correlations. Additive genetic correlations decreased steadily with increasing lag in age -yielding a rather smooth correlation surfaceto a minimum between birth and 280 days of 0.60 for PH and 0.70 for WOK. Previous, bivariate analyses considered birth and weaning weights, with average ages at weaning of 211 and 214 days for PH and WOK, respectively. Estimates of r A from an analysis fitting both genetic and environmental maternal effects (Model G2) were 0.69 for PH and 0.74 for WOK [26]. Corresponding estimates from the RR analysis at average ages agreed closely, with 0.67 and 0.71 for PH and WOK, respectively. Similarly, corresponding estimates of r P were 0.51 (PH) and 0.54 (WOK) for previous, bivariate analyses and 0.53 (PH) and 0.55 (WOK) for RR analyses.
Permanent environmental correlations. For most ages, estimates of r R decreased again with increasing time between records, correlations forming a plane descending steadily from the diagonal ridge at unity. For both breeds, however, there were "spikes" at the corners corresponding to the very youngest and very earliest ages. Presumably this was due to numerical problems -RR fitted for this random effect involved ages to the power 5. Higher order polynomials are notoriously problematic in this respect. Problems are exacerbated in areas involving fewest pairs of records and furthest from the mean. CFs have previously been observed to "misbehave" at the extremes of the ages considered [15,22,24], especially at higher -in particular unnecessarily highorders of polynomial fit.
Sizeable variances of higher order RR coefficients for permanent environmental effects due to animals (Tab. III) suggested though that the order of fit of k = 6 for this source of variation was not excessive. This was confirmed by a substantial decrease in log L when reducing the order of fit to a quartic polynomial (see Tab. II: k = 5253 vs. k = 5263 for PH and k = 5052 vs. k = 5062 for WOK). Whilst a stronger association between early post-natal weights which were least subject to maternal effects and weights at the latest ages, representing mostly post-weaning records, than with the intermediate weights is perceivable, estimates of r P did not show this pattern, indicating that the spikes were an artifact of combined effects of sampling variation, reduced variation at the highest ages, low numbers of records for pairs of extreme ages and high orders of polynomial fit.
Maternal correlations. The mean squared difference in estimates of permanent environmental maternal correlations (r C ) between breeds, calculated as above, was 0.0037. Estimates of r M for PH were more similar to estimates of r C in WOK than those of r C (in PH). When pooling genetic and environmental covariance estimates for PH and calculating an overall maternal correlation, the mean squared difference with r C for WOK was reduced to 0.00085. Estimates of maternal correlations formed a plateau close to unity from about two months of age (one month for r C for PH), indicating that maternal effects on weight from that age onwards were essentially identical. Estimates of r M and r C (for PH) between weights at birth and 211 days were 0.57 and 0.66, respectively, compared to values of 0.48 and 0.69 from previous analyses [26].

Additional analyses
As shown in Figures 1 and 2, calves tended to grow almost linearly up about 8 months of age with correspondingly increasing variances. Hence results from RR analyses which suggested that quartic or even quintic polynomials were necessary to model growth till just after weaning were surprising. Additional analyses for subsets of PH records were carried out to investigate possible causes for this phenomenon. Restricting ages in the data to 250 days and omitting any weights recorded after the actual weaning date (19 399 records), however, did not yield a reduction in the order of fit required. Scaled values for log L (+41 400) and BIC (−83 000) were −56. Eliminating birth weights and weights at very early ages (prior to 14 days of age) in addition, reduced the data to 16 438 records on 2 863 animals. However, it still did not allow the overall order of polynomial fit to be reduced. Scaled values for log L and BIC were 8.0 and 54.2 for k = 5263, −99.4 and 163.4 for k = 5243, −43.0 and 469.9 for k = 3263, and −191.2 and 260.7 for k = 3243, respectively, i.e. differences in information criteria between k = 5263 and models with lower orders of fit were somewhat reduced. Variances for the fourth and fifth RR coefficient for A and the fifth RR for R were considerably lower than in the complete data (see Tab. III), ranging from 0.5 to 0.8. An analysis for k = 5263 but omitting the fourth coefficient for A and the fifth coefficient for R reduced the number of parameters to 38, which decreased log L to −22.6 but reduced BIC to 9.9. Omitting the fifth coefficient for R for model k = 3263, however, did not prove advantageous (BIC 174.2), indicating that the quartic coefficient for A accommodated some variation due to the omitted quartic coefficient for R in the former analysis.
Estimates of σ 2 A for this subset were consistently lower than found for analyses including birth weight records, and up to 180 days of age agreed considerably better with univariate results. Estimates of the other variances were little changed, with a slight increase for the maternal components at the younger ages. This resulted in somewhat lower estimates of σ 2 P . Nevertheless, from 180 days onwards σ 2 P estimates were still markedly higher than for univariate analyses, similar to the pattern observed for WOK. As shown in Figure 6, this yielded markedly changed heritability estimates, in particular for the youngest ages. Similarly, estimates of m 2 and c 2 (not shown) were slightly increased, in particular for c 2 between 30 and 130 days of age estimates for the subset agreed very closely with corresponding univariate values. There was little difference in h 2 estimates for model k = 3263 and k = 5263, suggesting that a quadratic regression might suffice if birth weights were excluded from the RR analysis (or perhaps included in a bivariate analysis with birth weight as a trait with single records and a RR model for weights at other ages).
Yet, there appeared to be consistent, higher order environmental variation for direct effects. The data originated from a research station with a highly seasonal pattern of rains and thus pasture availability, winter rains being followed by almost complete summer droughts. As emphasized above, onset of summer conditions varied between years, with weaning dates varying accordingly from early November to early January, whilst calving dates remained constant. As noted for the analysis of mature weight records for cows in this experiment, analysis within contemporary groups removes systematic differences in weights, but not necessarily differences in seasonally induced variation [25]. It has to be assumed that the variation in higher order RR coefficients for direct, permanent environmental effects similarly reflects differences in seasonal variation.

Parameter estimates
Whilst there is a plethora of studies considering genetic parameters for birth and weaning weights in beef cattle (see, for example [17,29] for reviews), postnatal weights till weaning are seldom examined. Distinct differences between breeds in the importance of maternal effects have been reported. In particular, permanent environmental effects in Herefords are often found be almost twice as important than in other breeds, estimates of c 2 for weaning weight in the 20% to 30% range being common.
Results identify similar breed differences, with estimates of c 2 for PH above 20% between 40 and 230 days of age, while estimates for WOK were considerably lower. Breed differences in the importance of maternal environmental effects appear to be most important during the first 5 months. Moreover, highest values for c 2 occur earlier for PH (around 2-4 months of age) than for WOK (around 4-5 months). Large maternal effects in Herefords are commonly attributed to limited milk production in this breed. Means (± SD) for 14-hour milk production, measured by the weigh-suckle-weigh method in August each year, i.e. at an average age of calf of about 4 months, were 3.6 (±1.6) kg for PH and 4.9 (±1.9) kg for WOK [27]. It has been speculated that the substantially larger maternal effects in Herefords are due to an earlier decline of lactation curve than in other breeds. Differences between PH and WOK in age at which highest c 2 estimates were found are consistent with such hypothesis.
A series of univariate analyses of early weights in Brazilian Nelore cattle, based on large numbers of records, identified similar trends in estimates of genetic parameters. In particular, there was a corresponding decline in direct heritabilities after birth and a peak in the importance of maternal effects prior to weaning [2]. Estimates of m 2 varied little from 2 months of age till weaning. Ignoring a slight increase after 250 days, m 2 was highest at weaning, indicating that weaning weight recorded at about 200 days of age represented the most heritable selection criterion for maternal ability.
High order polynomial regressions were required to model changes in variability with age adequately, in particular for direct, permanent environmental effects. It was suspected that this could be due to the fact that a proportion of the latest weights represented post-weaning weights, but eliminating such records did not yield to a reduction in order of fit necessary. On the other hand, low order polynomial regressions are known to provide a poor fit for trajectories with a steep initial increase which then level off to approach an asymptote, such as growth curves ( [5], p. 102). Whilst the part of the growth curve considered here (birth to weaning) was still far from the asymptote represented by mature weight, this intrinsic behaviour of polynomial regressions may have contributed to the high order of fit needed.
Problems with overestimates of additive genetic variances were apparent for PH but not WOK when birth weights were included in the analysis. Both herds were managed in the same way and both data sets had a similar structure. Why such problems should arise for PH but not WOK is not clear.

Model selection
There has been concern that use of LRT to determine the "best" model to fit the data might favour overparametrised models. This lead to the use of information criteria -sometimes also referred to as penalised likelihoodswhich adjust for number of parameters estimated and sample size. As shown in Figures 3 and 4, estimates of variance components and genetic parameters from models selected on the basis of BIC differed only little from those obtained by models with the significantly highest log L (or minimum AIC) while reducing the number of parameters by 17 for PH and 13 for WOK. This suggests that BIC or related criteria should be employed rather than LRTs in selecting models for random regression analyses involving multiple random effects.

Alternative models
Results indicate that a relatively high order of fit is required to model individual variation in growth to just over 9 months of age through regression on orthogonal polynomials. Separation of direct and maternal, genetic and environmental effects required as many as three to four random effects to be considered. Fitting a quadratic to quintic polynomial for each and estimating the corresponding covariances among RR coefficients resulted in a large number of parameters to be estimated. Consequently, computational requirements for the analyses presented were extensive, most requiring weeks to complete. Furthermore, problems of sampling variation and correlation and in separating individual components were evident.
Whilst capable of modelling variation in growth adequately, polynomial regressions may not be the most appropriate functions. Standard growth curves, e.g. Gompertz functions or Brody's curve, however, are more suitable for data including mature weights and often don't fit well early in life. Future research should investigate alternative, more parsimonious models. On the one hand, these could include parametric curves (linearised if applicable) or splines, similar to applications to model test day yields in dairy cattle (see [4] and [42], respectively). On the other hand, there may be scope to reduce the number of parameters by modelling the within animal, permanent environmental covariance structure invoking a parametric correlation function described by one or two parameters in conjunction with a variance function to allow for heterogeneous variances (e.g. [6,32,33]).

CONCLUSIONS
Random regression models allow changes in traits and their variances over time to be modelled, and thus appear preferable for the genetic evaluation of beef cattle over the current multivariate approach which, somewhat arbitrarily, defines records to be different traits for different ranges in age at recording. However, this requires accurate estimates of covariance function for both direct and maternal genetic effects. This study presented a first attempt at estimating covariance functions for maternal effects on early growth of beef cattle.
Random regression analyses of early growth data, separating direct and maternal effects, are feasible albeit computationally very demanding. Furthermore, analyses were affected by numerical problems in fitting high order polynomials, and problems of estimating a large number of highly correlated parameters. As such, random regression analyses of growth data subject to maternal effects cannot yet be recommended as a routine procedure.
Regression on orthogonal polynomials proved well capable of modelling changes in the trait analysed and it's variability with time, although an unexpectedly high order of fit was necessary. Variance functions provided an effective and parsimonious way to accommodate measurement error variances increasing with age.
Results identified breed differences in the magnitude and peak time of importance for maternal effects, but very similar patterns of correlations. Implications for a RR model genetic evaluation scheme are that these might be modelled through common correlation functions whilst differences in variances could be accounted for through breed specific variance functions.

APPENDIX
Average information (AI) REML algorithms described in the literature generally consider the case where the variance matrix (V) of the vector of observations is linear in the parameters to be estimated. In that case, derivatives of V consist of sparse matrices with non-zero elements of unity and, more importantly, second derivatives of V are zero and the exact average of observed and expected information is readily calculated. If this is not the case, Gilmour et al. [9] suggest to approximate the exact average by a "simplified average" which approximates the second derivatives of the data part of the log likelihood with its expectation (which is the same asymptotically). Computationally, this is equivalent to ignoring the extra terms due to non-zero second derivatives of V.
Let R = Diag{σ 2 i } denote the diagonal matrix of measurement error variances pertaining to the vector of observations, with i = 1, . . . , N, and the i-th weight recorded at age a ij . As outlined in Meyer [20], terms required in an AI-REML algorithm are the first derivatives of R −1 and of log |R|. For R modelled by a VF as given in (4) and (5), with a * ij denoting the standardised age at recording for the i-th observation, these are: Variance Function To access this journal on line: www.edpsciences.org