Estimating genetic covariance functions assuming a parametric correlation structure for environmental effects

A random regression model for the analysis of "repeated" records in animal breeding is described which combines a random regression approach for additive genetic and other random effects with the assumption of a parametric correlation structure for within animal covariances. Both stationary and non-stationary correlation models involving a small number of parameters are considered. Heterogeneity in within animal variances is modelled through polynomial variance functions. Estimation of parameters describing the dispersion structure of such model by restricted maximum likelihood via an "average information" algorithm is outlined. An application to mature weight records of beef cow is given, and results are contrasted to those from analyses fitting sets of random regression coefficients for permanent environmental effects.


INTRODUCTION
Random regression (RR) models have become a preferred choice in the analysis of longitudinal data in animal breeding applications. Typical applications have been the analysis of test day records in dairy cattle and growth or feed intake records in pigs and beef cattle; see, for instance, Meyer [25] for references.
RR models are particularly useful when we are interested in differences between individuals, as we obtain a complete description of the trajectory, i.e. "growth curve", over the range of ages considered. A popular model involves regression on (orthogonal) polynomials of time. This does not require prior assumptions about the shape of the trajectory. Such RR have proven to be well capable of modelling changes in variation due to distinct events, such as weaning in beef cattle [25], or seasonal influences, e.g. [24]. However, frequently this required high orders of polynomial fit, and thus a large number of parameters to be estimated, accompanied by extensive computational requirements and numerical problems inherent to high order polynomials.
More generally in the analysis of longitudinal data, within-subject covariances between repeated records are often assumed to have a parametric correlation structure. In the simplest case, this might require a single parameter to specify correlations between records together with other parameters to model variances of records. A well-known example is the so-called "auto-correlation" structure. Other models involving a single parameter or low numbers of parameters (2,3,4) to model a correlation function are available, e.g. [2,11,29,31,40].
Pletcher and Geyer [33] presented an application of such models in the estimation of genetic covariance functions for age dependent traits in Drosophila. Their approach teamed a polynomial variance function (VF) to model changes in variances with age with a one-parameter correlation function (RF) to model correlations between different ages. However, estimation of the resulting covariance function (CF) used a reparameterisation of the covariance matrix among all ages in the data, as used by Meyer and Hill [27]. This resulted in computational requirements proportional to the number of ages in the data. Hence, their procedure is not readily applicable to large data sets arising from animal breeding applications with numerous different ages.
Recently, Foulley et al. [5] described an Expectation-Maximization type restricted maximum likelihood (REML) algorithm to estimate the covariance parameters for a model which combined a RR approach to model variation between subjects (e.g. genetic) with a single parameter RF to describe within subject covariances between repeated records. Their model included up to three parameters to model the latter, namely the parameter for the RF, the within subject variance and a measurement error variance.
Simple correlation models like those considered by Pletcher and Geyer [33] and Foulley et al. [5], generally imply stationarity, i.e. that the correlation between observations at any two times depends only on the difference between them, the "lag", not the times themselves. This might not be appropriate for animal breeding applications. Non-stationary correlation or covariance models are available, but usually involve more parameters. A common model, available in standard statistical analyses packages, is the so-called "ante-dependence" model [15]. A more parsimonious variant are structured ante-dependence models [32]. Pourahmadi [34] recently considered such models in a general mixed model framework.
This paper outlines REML estimation for RR models in animal breeding applications, assuming a parametric correlation structure for within animal covariances between repeated records. Both stationary and non-stationary models are considered. A numerical example comprising the analysis of mature weight records of beef cows is presented.

Random regression model
RR models commonly applied in animal breeding include at least two sets of RR coefficients for each animal, representing direct, additive genetic and permanent environmental effects, respectively. Let y ij denote the j-th record for animal i taken at time t ij . Assume we fit RR on orthogonal polynomials of time or age at recording. The RR model is then (1) with F ij denoting the fixed effects pertaining to y ij (often including a fixed regression on polynomials of time at recording), α im and γ im the additive genetic and permanent environmental RR coefficients for animal i, respectively, k A and k R the corresponding orders of polynomial fit, φ m (t ij ) the m-th orthogonal polynomial of time t ij (standardised if applicable), and ε ij the temporary environmental effect or "measurement error" affecting y ij .
Let α i = {α im } and γ i = {γ im } denote the vectors of RR coefficients for animal i of length k A and k R , respectively. Assume a multivariate normal distribution of records y ij , and mn } the matrices of covariances among RR coefficients, and σ 2 ε k the variances of measurement errors. Further, let y i be the ordered vector of observations for the i-th animal (ordered according to t ij ), and y of length M represent the complete vector of observations for all animals in the data, i = 1, . . . , N. Assume relationships between animals are known and taken into account, incrementing the number of animals in the analysis through inclusion of parents without records to N A . Let b of length N F denote the vector of fixed effects to be fitted with design matrix X, and α of length k A × N A and γ of length k R × N the vectors of additive genetic and permanent environmental RR coefficients. Design matrices for α and γ have non-zero elements φ m (t ij ), i.e. orthogonal polynomials evaluated for the times at which measurements are recorded.

K. Meyer
Let φ of size M×k R N denote the matrix of orthogonal polynomials evaluated for the ages in the data, with non-zero block of size n i × k R for the i-th animal. This is the design matrix for γ. The corresponding matrix for α is φ A of size M × k A N A , augmented by columns of zero elements for animals without records. Finally, let ε denote the vector of measurement errors corresponding to y. This gives Let y and γ be ordered according to animals, α be ordered according to RR coefficients. With A denoting the numerator relationship between animals and I N an identity matrix of size N, this gives With measurement errors assumed uncorrelated, Σ ε is diagonal and the mixed model equations and matrix (MMM) pertaining to (2) can be set up as for univariate analyses. Moreover, if Σ ε = σ 2 ε I M or Diag σ 2 ε d k , σ 2 ε can be factored from the MMM and be estimated directly from the residual sum of squares [22].
The covariance function due to permanent environmental effects of the animal (R) is estimated through K R . With R generally fitted to reduced order, i.e. k R smaller than the number of ages in the data, the resulting estimate of R, the permanent environmental covariance matrix among observations, is smoothed and has reduced rank. However, it does not have a pre-imposed structure. Whilst it is straightforward to estimate K R assuming a certain structure, this does not translate readily to R.

Equivalent model
We are, however, more interested in imposing a structure on R than K R . This can be achieved by fitting an equivalent model to (2) (4) with e of length M the vector of total environmental effects, i.e. the sum of permanent effects due to the animal and measurement errors. This has variance Like R, R * is blockdiagonal for animals. Permanent environmental covariances between records taken on the same animal are modelled through non-zero offdiagonal elements in the i-th block of R * , R * i . The MMM for (5) can be set up for one animal at a time, as for standard, non-RR multivariate analyses.
They can be thought of as derived from the MMM for (2) by absorbing γ, and computational requirements to factor the MMM are the same for both models.
Choosing (5) rather than (2), however, offers a much wider choice of parameterisation for R * and R, and allows for a chosen structure of R to be imposed easily.

Parametric correlation structures
Decompose R into the product of standard deviations and correlations R j the diagonal matrix of permanent environmental variances pertaining to y, and C = c j k the corresponding matrix of correlations. C is blockdiagonal for animals.

Variance function
Heterogeneous variances have been modelled through VF, e.g. [6,33,34], and this has been applied to measurement error variances in RR analyses [12,25,35]. Similarly, we can model the j-th element of Σ R or Σ 1/2 R as a function of the age at recording t ij . This can be a step function or, as more commonly used, a polynomial function, For instance, (9) with σ 2 R 0 the variance at the intercept, β r the coefficients of the VF and v the order of polynomial fit. Either variances (w = 2) or standard deviations (w = 1) can be modelled in this way. Functions (8) and (9) are advantageous when variances increase exponentially with time. In addition, they require less restrictions on the parameters of the VF than (7) to ensure that σ 2 R j > 0 for all j = 1, . . . , M. Whilst functions shown above involve ordinary polynomials (as in previous applications), use of orthogonal polynomials of t ij may be preferable to reduce sampling correlations between β r and thus improve convergence when estimating these parameters. Alternatively, for applications where variances show some periodicity, e.g. due to seasonal influences, a VF involving both polynomial and trigonometric terms [4] may be beneficial. In other instances, segmented polynomials [7] may be able to model changes in variances with time with fewer parameters.

Correlation function
Correlations between observations at different ages can be modelled as a function of the ages and one or more parameters of the RF. Correlation functions are stationary if a correlation between a pair of records depends only on the differences in ages at which they were taken -or lag-rather than the ages themselves. Most popular RF, including those given by (11) to (19) below, fall into this category.

Compound symmetry
In the simplest case, correlations between all observations for an animal (at different ages) are assumed to be the same.
with ρ a correlation, i.e. −1 < ρ < 1. This pattern is generally referred to as uniform correlation or compound symmetry (CS), and is the correlation structure assumed in the standard "repeatability model" analyses often used in the analysis of animal breeding data.

Auto-correlation
Let j k = |t ij − t ik | denote the lag in ages for a pair of records (y ij , y ik ) on the i-th animal. The so-called power, serial or auto-correlation function is then c j k = ρ j k (11) with −1 < ρ < 1 as above. This is the correlation structure generated by a continuous-time, first order auto-regressive (AR(1)) process.

Exponential model
An alternative way to model the correlation structure given by (11) is the exponential (EXP) model c j k = e −θ j k (12) with θ = − log(ρ) > 0. Again this parameterisation can be advantageous in terms of estimation, as it does not require the parameter of the RF to be constrained to an interval.

Gaussian model
In some instances, the decline in correlation with increasing lag is steeper than can be modelled with an exponential function of j k . In this case, the so-called Gaussian (GAU) exponential model which uses 2 j k may be more appropriate.
Diggle et al. [3] emphasize that in contrast to EXP, GAU is differentiable at j k = 0, and that for a sufficiently small time scale GAU has smoother appearance than EXP.

Other single parameter functions
Other RF involving different distributions but only a single parameter have been examined by Pletcher and Geyer [33]. All yield correlations which decrease with increasing lag. For instance, c j k = cosh(πθ j k /2) −1 (15) c j k = sin(θ j k )/ θ j k (16) c j k = 1 − cos(θ j k ) / θ 2 2 j k (17) are RFs based on the Cauchy distribution, the hyperbolic cosine, the characteristic function of the uniform distribution, and the characteristic function of the triangular distribution, respectively.
"Damped" exponential model A more flexible model can be obtained by adding a second parameter κ. This is a scale parameter which allows the exponential decay of the auto-correlation function to be accelerated or attenuated. Muñoz et al. [29] presented this for the serial correlation model c j k = ρ κ j k (18) pointing out that for κ = 1, κ = 0 and κ = ∞, (18) reduces to the serial correlation, compound symmetry and first-order moving average model, respectively. Alternatively, (12) can be expanded to [11]. Pletcher and Guyer [33] consider this as RF based on the characteristic function of the general stable distribution, with restriction 0 < κ < 2. In the following, (19) is referred to as damped exponential (DEX) model.

Other two parameter functions
A model which is not a special case of DEX, is the RF generated by a second-order auto-regressive process, which has parameters determined by the correlation between ages with lags 1 and 2 [29].
Any of the above RF ((11) to (19)) can be modified to allow for a proportion τ of the correlation independent of age effects [11] (20) with c * jk a function of the lag in ages as modelled above, and τ estimated as an additional parameter (yielding a three-parameter RF if extending (19)).

Structured ante-dependence model
Another class of models employed in the analysis of "repeated" records or longitudinal data are the so-called ante-dependence (AD) models, e.g. [15]. These are loosely related to time series models, in such that the j-th record on an animal depends on and is correlated to a number of its predecessors [3]. In contrast to the parametric correlation structures considered so far, AD models allow for non-stationary correlations.
For an AD model of order s, AD(s), a record y ij in the ordered vector of observations y i for animal i is assumed to depend at most on records y i ( j−1) , . . . , y i ( j−s) , but to be independent of any other preceding observations y i ( j−s−1) , . . . , y i 1 . This yields a correlation matrix with elements on the first s subdiagonals as variables, and the elements of the remaining subdiagonals (s + 1, . . . , n − 1) determined by the former. Consequently, the corresponding inverse is a banded matrix, with only the elements of the leading diagonal and first s subdiagonals being non-zero [15]. Hence, for n different times of recording, an unstructured AD(s) model has (s + 1)(2n − s)/2 parameters, n variances and sn − s(s + 1)/2 correlations.
For s = 1, a first-order AD model, there are n − 1 correlations on the first sub-diagonal of the correlation matrix, c j ( j+1) . The other correlations are given by a simple multiplicative relation for j = 1, n and k = j + 2, n (21) [31]. For s > 1, the functional relationship with the elements of the first s subdiagonals is more complicated. In that case, a parameterisation in terms of the inverse of the corresponding covariance matrix -also called the "concentration matrix" of the AD -or it's Cholesky decomposition is often preferred. Whilst an AD(s) model with low s has considerably less parameters than a full multivariate, unstructured model (which has n(n + 1)/2 parameters), it can still involve impractically many parameters. Structured ante-dependence (SAD) models [32] assume a functional relationship between the parameters of an AD model, and thus provide a more parsimonious representation. Firstly, variances are considered to be a function of the time at measurement, with the function involving a small number of parameters. Zimmerman et al. [39], Núñez-Antón and Zimmerman [32] and Pourahmadi [34] consider polynomial VFs as described above (see (7) to (9)). Secondly, the correlations on the first s subdiagonals are determined by the times of recording and 2s parameters, ρ k and κ k , respectively: for k = 1, s and j = k + 1, n (22) [32] with 0 < ρ < 1 and for κ k = 0 (23) Function (23) applies a deformation (Box-Cox power transformation) to the time scale which facilitates non-stationarity of correlations. For κ < 1 equidistant correlations are increasing with age. Conversely, κ > 1 implies lower correlations between records with equal lag at higher ages [39]. For s = 1 and κ = 1, the RF (22) reduces to (11), the auto-correlation function.

ESTIMATION OF COVARIANCE AND CORRELATION FUNCTIONS
Parameters of covariance, correlation and variance functions are readily estimated by restricted maximum likelihood (REML). This may involve a derivative-free procedure, an "Average Information" (AI-REML) algorithm [8] or an Expectation-Maximization (EM) algorithm, as described by Foulley et al. [5]. Various authors consider REML estimation in the analysis of longitudinal or spatial data, but often do not go further than specifying the log likelihood and using a simple search procedure, such as the simplex method of Nelder and Mead [30], to locate its maximum, e.g. [3,31,39]. Others describe maximum likelihood estimation using Newton-Raphson type algorithms, e.g. [13,18,29]. Gilmour et al. [8] consider AI-REML estimation for models with correlated residuals in a general formulation.

The likelihood
The REML log likelihood for (4) is −2 log L = const + log |G| + log R * + log |C M | + y Py (24) where C M is the coefficient matrix in the mixed model equation pertaining to (4) and y Py is the sum of squares of residuals. Both y Py and log |C M | can be evaluated simultaneously as described by Graser et al. [9], by factoring the corresponding MMM M is large but sparse, with N M = N F + k A N A + 1 rows and columns. For R * blockdiagonal it can be set up for one animal at a time, as for corresponding multivariate analyses. Factoring M into LL with L a lower triangular matrix with elements l ij (l ij = 0 for j > i) gives log l kk and (26) y Py = l 2 e.g. [28]. The other components of (24) can be evaluated as This involves determinants of small matrices only, of size k A and the number of records for each animal, respectively. For some correlation structures, closed forms for the corresponding inverse correlation or covariance matrices and determinants exist. In some cases, in particular for analyses assuming Σ ε = 0, this can be exploited to reduce computational requirements to evaluate (29).

AI-REML algorithm
Maximisation of log L via AI-REML requires first derivatives of (24) and the average of observed and expected information [8]. The latter is proportional to second derivatives of the data part, y Py, of the likelihood. These can be determined as for standard multivariate analyses, using sparse matrix inversion and repeated solution of the mixed model equations [19] or automatic differentiation of the MMM [20].

First derivatives
Derivatives of log |C M | and y Py can be determined through automatic differentiation of the Cholesky factor of M, as described by Smith [36]. This requires the derivatives of M with respect to the parameters to be estimated.

∂M ∂K
for covariances among the RR coefficients for additive genetic effects, α i . The derivative ∂K A /∂K A mn has elements of unity in position mn and nm and zero otherwise. Hence, with K mn A the mn-th element of K −1 A and δ mn Kronecker's delta, i.e. δ mn = 1 for m = n and 0 otherwise. Similarly, and θ R standing in turn for the parameters of the permanent environmental VF (P ), σ 2 R0 and regression coefficient β r , the parameters of the CF, ρ or θ and κ, and the parameters of the VF for measurement error variances (E ).
Let ∂L/∂θ R , with elements ∂l ij /∂θ R , denote the derivative of L with respect to θ R obtained by differentiation of L as described by Smith [36]. First derivatives of the last two terms in (24) are then [28] ∂y Py

K. Meyer
Derivatives of the other two terms in (24) can be evaluated indirectly For R * linear in the parameters to be estimated, (36) simplifies analogously to (35), see [28].

Second derivatives
AI-REML algorithms [8,14,19,20] have generally considered the case where V is linear in the parameters to be estimated, as, for instance, for standard multivariate analyses. This gives second derivatives of V which are zero, and the average of "observed" and "expected" information for parameters θ r and θ s is equal to − 1 2 y P ∂V ∂θ r P ∂V ∂θ s Py (37) For cases where V = p θ p (∂V/∂θ p ), as for our parameterisation of the residual covariance matrix in terms of a variance and parametric correlation function, second derivatives of V are non-zero. For such models, Gilmour et al. [8] suggest to approximate the exact average by a "simplified average" information. This is derived by approximating ∂ 2 y Py/∂θ r ∂θ s by its expectation. Asymptotically the two are the same. Computationally, this is equivalent to ignoring extra terms involving non-zero second derivatives of V.
Rewrite (37) as b r Pb s with b r = ∂V/∂θ r Py. For θ r = K A mn , (38) with φ A q and α q denoting the sub-matrix of φ A and subvector ofα, respectively, for the q-th RR coefficient. Similarly, for θ r a parameter of where " + " denotes the direct matrix sum. As shown previously [20], crossproducts b r Pb s can be evaluated by replacing y in (25) with B = b 1 |b 2 | . . . |b p , a matrix containing one column b i for each parameter, and expanding the MMM to (40) Factoring M B then overwrites B PB with elements b r Pb s . With the Cholesky factor of C M already evaluated (in factoring M), this is computationally undemanding.

Derivatives of R *
Evaluation of Q R (33), the first derivatives of R * (36) and vectors b r pertaining to parameters of R * (39) requires the partial derivatives of R * with respect to the parameters to be estimated, as well as products and traces involving the inverse of R * . Corresponding terms for the parameters of a polynomial VF for measurement error variances under model (2), i.e. the simple case of a diagonal residual covariance matrix, have been given by Meyer [25].

Correlation function
Derivatives of c j k need to be determined for each RF separately. For instance, for the auto-correlation function, (11) or (18) Similarly, for an exponential model, EXP (12), GAU (13) or DEX (19), with fixed values for κ where appropriate.
(a) w = 1 to model standard deviations, w = 2 to model variances.

Variance functions
Derivatives of the diagonal matrix of measurement error variances are with θ ε standing for the parameters modelling changes in temporary environmental variances over time. Expressions for parameters of a VF, E , are analogous to those given below for permanent environmental variances. For the parameters of P , θ R = σ 2 R0 or σ R0 and β r for r = 1, . . . , v, Derivatives of σ R j depend on the function and parameterisation chosen. Values of ∂σ R j /∂θ R for functions (7) to (9) to model either standard deviations or variances are summarised in Table I.

Data
Data consisted of January weights of Polled Hereford beef cows between 2 and 10 years of age, analysed previously fitting RR on Legendre polynomials of age for both additive genetic and permanent environmental effects of the animal [23].
Records originated from a selection experiment carried out at the Wokalup research station in the South West of Western Australia; see [26] for details. With a Mediterranean climate, characterised by Winter and Spring rains and pasture growth, and subsequent almost complete drought in Summer and Autumn, January weighings tended to record cows at their top weight during the year. Short mating periods resulted in the bulk of calves born in April and May each year. Ages at weighing for records selected thus ranged from 19 to 119 months, with up to 9 weights per cow. In total there were 3 320 records on 850 cows, offspring of 149 sires and 423 dams, with 165, 153, 103, 124, 83, 75, 47, 50 and 50 animals having 1, 2, . . . , 8 and 9 records available.

Analyses
Data were analysed assuming a parametric correlation structure for covariances between records taken on the same animal. Models CS, EXP, GAU, DEX and SAD were considered, teamed with polynomial functions to model permanent environmental standard deviations ( (7) with w = 1) of order v = 0 to 7. Measurement error variances were in turn considered homogeneous (e = 1), to be represented by a polynomial function ( (7) with w = 1) as for permanent environmental effects, or to change with year of age, fitting a step function with e = 7 classes (2, 3, 4, 5, 6, 7, 8-10 years) as in previous analyses [23]. In the following "RF.VFv.Ee" denotes an analysis fitting correlation function RF with polynomial function P for permanent environmental standard deviations to order v, and e measurement error variances. A model fitting a polynomial function E for temporary environmental standard deviations with v coefficients is described as "RF.VFv.E1+VFv ". For example, "SAD.VF3.E1" represents an analysis assuming a structured ante-dependence model for the correlations between records for an animal with a cubic variance function to model changes in permanent environmental variances over time, and temporary environmental effects which are considered to have homogeneous variances.
For comparison, data were also reanalysed assuming permanent environmental covariances followed a pattern as described by fitting a set of RR coefficients on Legendre polynomials (LP) of age at recording, i.e. R = φK R φ . As in previous analyses, estimates of K R were forced to have reduced rank, by setting eigenvalues less than 0.001 to zero, thus reducing the number of parameters to be estimated accordingly. In the following, "LP.Rk R rk R .Ee(+VFv )" denotes an analysis fitting LPs for permanent environmental effects to order k R , with estimated covariance matrix of rank k R (other terms as above).
Analyses were first carried out on a "phenotypic" level, assuming all covariances between records were due to animals'permanent environmental effects. This facilitated computationally undemanding examination of a wide range of variance functions and correlation structures. Secondly, analyses allowed for (co)variances between individuals by fitting RR coefficients on LPs of age K. Meyer due to animals' additive genetic effects, incorporating all pedigree information available. "Ak A rk A " denotes an analysis fitting LPs to order k A with estimated covariance matrix K A of rank k A .
As in previous analyses, fixed effects fitted included contemporary groups, defined as year-paddock subclasses, and a cubic regression on Legendre polynomials of age. Parameters of variance, correlation and covariance functions were estimated by REML, using a combination of derivative-free and AI-REML algorithms. Analyses were carried out using program DXMRR 1 [21], extended to accommodate model (4) with parameterisation of R * as described above.
Results from different analyses were compared by examining estimated variances and correlations for ages represented in the data. In addition, maximum values of log L and the REML forms of Akaike's (AIC) and Schwarz' (BIC) information criterion, e.g. [38] were contrasted.

"Phenotypic" analyses
Results from analyses ignoring variation between animals are summarised in Table II. Values of log L are clearly dominated by the order of fit for VF P and the degree of heterogeneity allowed for measurement error variances. For e = 1, P needed to be a cubic or higher order polynomial, for estimates of variances at mean age not to be drastic over-or underestimates. Estimated correlations between records one month apart were close to unity. At equal order of fit for P , v, two-parameter RFs (DEX and SAD) yielded higher log L than single parameter RFs (CS, EXP and GAU), but there was no advantage of the non-stationary SAD over the stationary damped exponential correlation structure.
Assuming homogeneous σ 2 ε , parametric RF had higher log L than analyses fitting LP involving similar numbers of parameters, presumably because more parameters were available to model changes in variation with time. For e = 7, however, there was little difference, with the heterogeneous σ 2 ε accounting for differences in variation not modelled by covariance function R. Model "DEX.VF3.E7" with 13 parameters had minimum AIC, though the corresponding log L was not very different from that for "LP.R4r2.E7". Involving a more stringent penalty for the number of parameters fitted, BIC suggested that "EXP.VF3.E1" with only six parameters sufficed to model the covariance structure adequately.
Estimated variances for selected analyses are shown in Figure 1, and corresponding correlations are given in Figure 2. For e = 1 and v ≥ 4 or v = 3 together with heterogeneous σ 2 ε , estimates of the total variance differed little between models of analyses. Differences were largest at the highest ages which were represented by comparatively few records only. For DEX, however, estimates of σ 2 ε were consistently lower than for the other models, in particular if a VF E was fitted or e > 1.
As shown in Figure 2, correlation functions EXP and DEX force estimates of within-animal correlations between repeated records to decrease steadily with increasing lag. Corresponding phenotypic correlation (r P ) estimates follow a similar pattern, but fluctuate if variances are heterogeneous.    Hence average values together with their ranges are shown for r P . Whilst EXP and DEX imply stationarity for r R , analyses fitting LPs do not impose this restriction, and values for r R shown are averages (with ranges), as for r P . Previous analyses fitting LPs [23] had reported an inexplicable peak in estimates or r R between the youngest and oldest ages, and attributed this to sampling variation. This peak is reflected in increases of the average estimate of r R for lags of more than 70 months for analyses "LP.R4r3.E1" and "LP.R4r3.E7". Similarity in values of log L and estimates of variances for the latter and "DEX.VF3.E7" suggests that differences in estimates of r R are indeed spurious.  Table III gives results from analyses fitting a set of RR coefficients for animals' additive genetic effects in addition. Clearly, there is variation between individuals which should not be ignored. Likelihoods increased up to k A = 4, the order of fit identified previously [23]. Fitting K A with rank k A = 2 was sufficient throughout. Again, BIC favoured a model with homogeneous σ 2 ε , "A3r1.DEX.VF5.E1" with 12 parameters, while AIC selected the model of previous analyses, "A4r2.LP.R4r2.E7" with 21 parameters.

Genetic analyses
Estimates of variances and correlations for the "best" models and others closely related to them are shown in Figures 3 and 4, respectively. For comparison, corresponding estimates from a "standard" multivariate analysis treating observations at each year of age (2 to 10) as a different trait are given in Figure 5. Most notable again is the similarity in estimates of the total variance for all analyses shown while the partitioning into invidual sources of variation differs. Problems with sharply diverging estimates of genetic and permanent environmental variances for the highest ages in this data set were observed before and attributed to the influence of atypical records at the upper extreme.
Similarly, estimates of the average r P (Fig. 4) for the analyses shown agree well, ranging from 0.7 to 0.8 to about 0.4. Estimating K A with rank k A = 1 forced all estimates of the genetic correlation (r A ) between records at different ages to be unity, while k A = 2 allowed for a decline in r A with increasing lag.   No. of parameters.

DISCUSSION
Whilst analyses of longitudinal, spatial or similar data assuming a parametric correlation structure or covariance function are commonplace in other areas of applied statistics, they have found few applications in the analysis of animal breeding data, e.g. [1,37]. Harville [10] considered auto-regressive random effects in linear mixed models, and suggested that modelling of dairy records might be improved by considering permanent environmental effects of cows as auto-regressive rather than constant across lactations. "Random coefficient" models are often found to be of little advantage and somewhat cumbersome, the resulting covariance functions generally yielding less readily interpretable results than the parameters of a parametric correlation function. Moreover, low degree polynomial regressions tend to provide a poor fit for trajectories with steep initial increases which then level off to gradually approach an asymptote [3], p. 102, i.e. high order polynomials and corresponding numbers of parameters are often required to model the covariance structure for growth curve type of analyses.
In contrast, animal breeders have embraced random regression models for the analysis of longitudinal data, in particular test-day records of dairy cows and growth data for pigs and beef cattle. This can be attributed to several factors.
Firstly, quantitative genetic analyses are invariably concerned with the variation between animals, while other areas of statistics are often content with modelling within-subject covariances only. Fitting a set of additive genetic RR coefficients provides estimates of genetic merit for the whole range of ages considered, and allows ranking of animals to change with time. RR models are thus an obvious choice if we are concerned with (genetic) differences between individuals. Assuming a RR model and resulting covariance structure on a genetic level, it seems natural to apply the same model to other random effects, for instance permanent environmental effects due to the individual.
Secondly, parametric models are usually formulated assuming homogeneous variances. Variances of "repeated" record traits of interest to animal breeders, however, often change with time, if only due to scale effects. While extensions to heterogeneous models by replacing single variances with a (polynomial) VF are straightforward, these are seldom described. RR models account for changes in variances with time, and for RR involving orthogonal polynomials, do not require any specific assumptions about the shape of the resulting VF.
Thirdly, estimates of genetic covariance matrices arising from RR model analyses can be thought of as smoothed versions of corresponding estimates from an unstructured, multivariate analysis treating records at different ages as different traits. Covariance functions which give covariances between records at individual ages as function of orthogonal polynomials of the ages and the elements of a matrix of coefficients (K), have been described as "infinitedimensional" equivalent to covariance matrices in standard multivariate analyses. Estimates of the eigenvalues and eigenfunctions of such CFs can be obtained directly as the eigenvalues of K and from the corresponding eigenvectors. For genetic covariance functions, these statistics provide valuable insight into the effects of selection for the trait considered. See Kirkpatrick et al. [16,17] for further details.
Last, but not least, RR models provide a computationally feasible way to estimate CF for large data sets with records coming in at "all ages", as are typical for data from livestock recording schemes. Estimating K as the matrix of covariances between RR coefficients requires mixed model equations of size proportional to the number of regression coefficients to be manipulated, rather than proportional to the number of ages or even the number of records. In contrast, general algorithms for the analysis of longitudinal data fitting parametric CF or RF are frequently formulated in terms of the variance matrix of the vector of observations, i.e. involve the latter, e.g. [3], Chap. 4. This paper shows how the preferred RR approach for genetic effects can be combined with a parametric correlation model for within animal covariances. It extends the model suggested by Foulley et al. [5] to a variety of correlation structures, both stationary and non-stationary, and heterogeneous variances. An average information algorithm for the estimation of the parameters for the covariance, variance and correlation functions describing the dispersion structure for such model, by REML is outlined. It involves computations proportional to the number of animals and the order of polynomial fit for genetic effects. Extensions to models involving additional sets of RR coefficients for other random effects, e.g. maternal effects, or more complicated correlation functions are straightforward.
Application to a data set of mature weight of beef cows showed that assuming a parametric correlation structure for within animal covariances can result in a considerably more parsimonious model than a RR model for permanent environmental effects. The example showed further that parametric correlation functions can eliminate erratic and inexplicable estimates of correlations between records at extreme ages which have been encountered with a RR. Some differential partitioning of the total variance was evident though. This occurred in particular, when both temporary and permanent environmental variances were assumed heterogeneous. An alternative for this scenario might be the use of a link function as suggested by Foulley et al. [6], which would assume a functional relationship between the two sources of variation and thus might not only reduce the number of parameters required further but also alleviate problems of erratic partitioning. Further work is required to examine the effect of different assumptions about the structure of within-animal covariances on the partitioning of the phenotypic variance and, thus, their potential impact on estimates of genetic parameters and predictions of breeding values.