Computing approximate standard errors for genetic parameters derived from random regression models fitted by average information REML

Approximate standard errors (ASE) of variance components for random regression coefficients are calculated from the average information matrix obtained in a residual maximum likelihood procedure. Linear combinations of those coefficients define variance components for the additive genetic variance at given points of the trajectory. Therefore, ASE of these components and heritabilities derived from them can be calculated. In our example, the ASE were larger near the ends of the trajectory.


INTRODUCTION
Random regression (RR) has been widely used for genetic analysis of longitudinal data from many of the major animal breeding industries world wide and has also been implemented into routine large scale animal breeding applications [4]. Estimates of derived genetic parameters such as heritability at given points along the trajectory are commonly published from such studies and comment is often made about the accuracy and robustness of RR models. However, there have been no attempts to quantify the accuracy of such estimates for different parts of the trajectory from RR analyses using residual maximum likelihood (REML) methods. In contrast, Meyer [9] published confidence intervals of genetic parameter estimates derived from Bayesian analyses using Gibbs sampling. With REML estimation by the average information algorithm, approximate variances of variance components are obtained from the inverse of the information matrix. Variance components as well as heritabilities at given trajectory points can be calculated from variances of random regression coefficients and therefore approximate standard errors (ASE) of these derived parameters can also be obtained. The aims of this note are to describe how to calculate ASE for genetic parameter estimates derived from RR models and to apply the method to a field data set.

Random regression model
Consider a variance-covariance (VCV) matrix G 0 of rank t for repeated measurements of weight at t given trajectory points (e.g. ages). Under the covariance function (CF) approach defined by Kirkpatrick et al. [5], G 0 is modelled with a reduced number of parameters. The genetic CF of order k, where k ≤ t, can be estimated from G 0 such that: whereĜ is an approximation of G 0 . Meyer [8] showed that K can be estimated directly from data using RR. The matrix K of order k contains the variance components for the RR coefficients in the model. The matrix Φ of order t × k contains orthogonal polynomial coefficients evaluated at t standardised trajectory points (ages) with elements φ ij = φ j (x i ), being the jth polynomial coefficient for the ith point x i [6]. The covariance structure for the environmental effects is fitted as an unstructured t × t covariance matrix. This yields the model: where y i is the vector of t i observations measured on animal i, b is a vector of fixed effects, α i a vector of additive genetic RR coefficients and e i a vector of residual errors pertaining to y i . X i and Z i are design matrices relating b and α i to y i , where Z i contains the elements Φ pertaining to ages in the data. Extending the model to n individuals, the corresponding variances are defined as var(α) = K ⊗ A, where K contains the additive genetic variances and covariances for the RR coefficients, A is the numerator relationship matrix among individuals and the symbol ⊗ denotes direct product. The solution for K can be used as in equation (1) to calculate the variances and covariances among defined trajectory points.

Calculation of standard error of parameters derived from RR coefficients
Consider a genetic variance covariance matrix,Ĝ, derived from equa- We can write the elements ofĜ in vector form, such that the variances and covariances of these parameters can be summarized in a matrix. Hence, equation (1) can also be written as is the vector form ofK of dimensions (k × k) × 1 achieved by stacking the columns ofK below one another, and similarly vec(Ĝ) is the vector form ofĜ of dimensions (t × t) × 1. It can be checked for a small example that equations (1) and (3) are equivalent, written in matrix and vector form respectively. The variance of estimates inĜ can be calculated in a similar manner whereby where var(vec(K)) has dimensions (k ×k)×(k ×k) and var(vec(Ĝ)) is a (t ×t) by (t × t) matrix. Var(vec(K)) can be approximated from the appropriate elements of the inverse of the average information matrix in a REML procedure (e.g. as given in the *.vvp file in ASReml) [3]. The same principles apply to other random effects in the RR model, and the covariance between variances of different random effects. Hence, this methodology can be extended to the matrices estimated for other random regression effects and the covariance between random effects. Subsequently these matrices are summed as in equation (5) to build a matrix containing estimates of variance of phenotypic (co)variance components, var(vec(P)), which also has dimension (t × t) × (t × t).
For functions of variance components (such as heritabilities) a Taylor series expansion can be used to approximate the variance of a variance ratio as detailed by Lynch and Walsh [7]. For the ratio of genetic to phenotypic variance, we get: whereĝ i,i andp i,i are elements ofĜ andP, var(ĝ i,i ), var(p i,i ) and cov(ĝ i,i ,p i,i ) represent variance and covariance of genetic and phenotypic variance at time i. The ASE for the heritability estimate at time i (for univariate and RR estimates) is obtained by taking the square root of equation (6).  A VCV matrix for additive genetic and phenotypic effects for weight over a 450-day trajectory was derived based on the analysis performed by Fischer et al. [2]. Data for this analysis originated from the LAMBPLAN database and consisted of 16 826 weight records on 5 420 Poll Dorset sheep. The number of records at different ages is represented in Figure 1.

Example of application of method to RR coefficients estimated from field data
Fischer et al. [2] used RR to estimate CF coefficients for direct and maternal genetic and environmental effects. The model also included heterogeneous residual variance across ages of measurement. ASReml [3] was used for this analysis. Based on a third order CF for additive genetic effects, a VCV matrix (Ĝ) was constructed for weights at 10 equidistant ages (i.e. defining Φ). Similarly, VCV matrices were derived for the other random effects. Furthermore, adding the resultant variance matrices together resulted in a phenotypic VCV matrix (P) with (co)variance components for weights at the 10 equidistant ages.
We then obtained the variance of vec(Ĝ) as in equation (4) and similarly for the two types of maternal effects, which in this case were all matrices of dimensions 100 × 100. Following equations (4), (5) and (6) we obtained the ASE of the heritability estimate for each age. Results from this example are shown in Figure 2.
In addition, a series of piecewise estimates at specific ages were obtained using the equivalent univariate model (direct and maternal genetic effects only, no permanent environmental effects fitted). Estimates were taken from day 50 onwards using a 50-day age window for each univariate estimate up to 500 days and these are also shown in Figure 2. As can be seen in Figure 2, the ASE for heritability estimates from RR are lower (0.05−0.07) in the middle of the trajectory, and higher (0.08−0.13) at the ends of the trajectory. The pattern at the end of the trajectory is largely due the nature of polynomials, which have no constraint at the ends. This result is consistent with that of Fischer and van der Werf [1] who demonstrated problems associated with polynomials, in particular high order Legendre polynomials.

RESULTS AND DISCUSSION
The ASE for heritabilities at given ages obtained using RR were smaller (0.04−0.13) throughout much of the trajectory than those obtained from piecewise univariate analyses of portions of the same data within 50 day age windows throughout the 450 day trajectory (0.04−0.15). The large standard errors for specific univariate estimates were largely a function of fewer data within that age window, in particular beyond 200 days. Furthermore, the larger ASE values at the ends of the trajectory are largely due to polynomial instability, which is further exacerbated by fewer data at these ages. Instability has been shown in heritability estimates at the ends of the trajectory when using polynomials in RR, even in cases where there was more data at the ends [1].
Furthermore, the ASE obtained using the RR model were comparable or lower than those obtained from piecewise univariate analysis except at the ends of the trajectory, showing where RR estimates are more accurate. It is particularly useful to obtain standard errors of heritability at the ends of the trajectory as concerns are often raised about the robustness of such estimates due to the polynomial nature of the covariance function and justification of this concern was demonstrated in this study. Moreover, the same method can be applied to other random effects to get standard errors for their respective proportion of phenotypic variance (e.g. maternal heritability).

CONCLUSION
This study demonstrated a method for computing ASE for genetic parameter estimates derived using RR models applied to a field data set. The method produced plausible standard error values for estimates of heritability and this provides insight into the discussion of robustness and accuracy of RR estimates of heritability at specific age points.