Estimating genetic parameters using an animal model with imaginary effects

Summary - The estimation of additive genetic variance by maximum likelihood is discussed. The extension of reduced animal models, when parents are not inbred to allow the use of existing algorithms for maximum likelihood estimation of additive genetic variance, is described. This involves the introduction of imaginary effects with negative variance, and leads to computation using complex arithmetic. Methods are developed to allow the computation to be carried out using only real arithmetic. This method has computational advantages when only a small proportion of animals have offspring. genetic parameter / animal model / estimation / maximum likelihood


INTRODUCTION
Additive genetic variance and heritability have most commonly been estimated in animal breeding data information from collateral relatives, such as half-sibs or fullsibs, or non-collateral relatives, such as parent-offspring. Covariances generated by these relationships provide the most information for estimating additive genetic variance. However, there is interest in combining these alternative estimates (Nicholas and Hill, 1974), and also using the additional information from other genetic relationships.
Linear models can be formed that contain genetic and environmental effects for each animal, commonly called individual animal models (eg, Quaas and Pollak, 1980). In theory at least, standard or restricted maximum likelihood estimation (ML or REML) procedures can be applied to individual animal models for variance component estimation (Harville, 1977). One problem is that this model generates as many equations as there are animals in the pedigree and so there have been attempts to reduce the computational effort. Thompson (1977) considered the case of two generations and showed how to develop REML estimating equations just for animals in the first generation. More generally, this estimation can usually be interpreted using predicted breeding values. Quaas and Pollak (1980) showed that a reduced animal model (RAM) has advantages in calculating predicted breeding values. Notably, only equations for animals with offspring are needed. Henderson (1986) and Sorensen and Kennedy (1986) showed that REML and minimum norm quadratic (MINQUE) estimation can be expressed with advantage using a RAM.
However, these (REML) methods are iterative and, as presented, require the inversion of a matrix with as many rows and columns as the number of parents in each round of iteration. For some variance component problems, this computational burden can be reduced by using orthogonal transformations to give simpler equations, either in terms of a diagonal matrix (Patterson and Thompson, 1971;Dempster et al, 1984), or in terms of a tridiagonal matrix (Smith and Graser, 1986).
These results are not directly applicable to estimation methods based on a RAM. This paper shows how a RAM can be modified by the introduction of extra random effects with negative variances, so that the resulting matrices involved in estimation can be tridiagonalized using existing algorithms based on Householder transformations. Complex numbers are needed to take account of negative variances, but a new algorithm is given that takes account of the special structure of the matrices involved, and reduces the need for arithmetic based on complex numbers.

THE MODEL
If additive genetic covariances are the only source of covariances between records, then a linear model for the individual animal model (IAM) can be written as: with E(y) = Xb, E(a) = E(e) = 0 and V(a) = Ao, 2, V(e) = Ia5 and cov(a,e') = 0, where y, b, a and e denote the vectors of observations, fixed effects, animal effects and residual errors. X and Z are the corresponding incidence matrices, and for simplicity, we assume X had r columns and rank r. With single records per animal, Z = I nt , where m denotes the number of observations. A is the numerator relationship matrix (NRM) between animals. As described by Quaas and Pollak (1980), with a RAM, the vector of animals is divided into parents, ie, animals which have offspring (subscript P in the following), and non-parents, ie, animals without progeny (denoted by subscript N). The additive genetic value for non-parents is then partitioned into contributions from parents and parts due to the Mendelian sampling. Including the latter, together with the residual e in eqn(l) into a new residual error, gives the RAM as a reparameterisation of eqn(1): Let p denote the number of parents and q the number of non-parents, then (p + q = m). Z N is then a matrix of order qxp, with elements z2! = 0.5, if j is a parent of i, and zero otherwise. This representation, strictly speaking, assumes all parents of non-parents are in the data and one record per individual. However, eqn(2) can be easily modified, if that is not the case. The residual error for a non-parent has variance at v = a5 + 0.5 QA , if both parents are known, and variance at v = QE +0.75 QA , if only one parent is known, provided they are not inbred. In the following, it is assumed that there is equal parental information for all non-parents and that parents are non inbred ie, that a 2 w is constant. Estimation equations, based on the RAM, usually have to be reconstructed after each variance iteration because the ratio of error to additive variance changes. Therefore, for computational reasons, it is desirable to have a model where the vector of residuals has homogeneous variance. This can be achieved by a further reparameterisation, adding effects to either parents or non-parents. As there are usually fewer parents than non-parents, fewer equations will be generated if eqn (2) is expanded to: The variance of e, is such that var(e j ) + var(e P -e j ) = var(ep) = Ipas and var(ep &mdash; e l ) = lo,2&dquo; so that ep &mdash; e j and e N have the same variance. Hence, var(e 1 ) _ -C 2 UA Ip, with c 2 = 0.75 and 0.50 for one and both parents known, respectively. Because of the assumption that there is equal parental information for all non-parents, c 2 is the same for all e, values. Normally, variances are assumed to be positive so there is a slight difficulty in interpreting e l with a negative variance.
However, if a D are effects with variance a A 2 Ip, then defining e j = ica D , where i = i , then the variance of e j = -C 2 0 , 2 A Ip. Hence, e i can be thought of as imaginary effects as they are a multiple of i.
Hence, eqn (3) (7) gives predictors of the effects ap and a D . Note that ap is real and a D is imaginary, so that ia R = a D with a R real, and the predictor for e l = -ca R is a real quantity, as would be anticipated. REML equations to estimate QA and Q yy are (in terms of real quantities): where r is the rank of (X' p Xp + X%XN ) , and n is the number of elements in both ap and a D .
These equations involve a 2 A and Q yy through A = !yy/QA on both sides of eqns (8) and (10), and have to be solved iteratively. To simplify eqns (8) and (10) and avoid the inversion associated with Ap l , Smith and Graser (1986) and Meyer (1987) suggest writing Ap as LL' and using a L = L-l ap in a model equivalent to eqn(3). The effects a L are uncorrelated. Using Meyer's results and eqns(8) and (10), gives an alternative form of the EM algorithm.
(A I is not to be confused with A or Ap, numerator relationship matrices). This iterative scheme requires the inversion of (A I + IA) in each iteration. Smith and Graser (1986) have shown that these computations can be reduced by writing A l as (PA!P') where P is an orthogonal matrix and A* is a tridiagonal matrix. Using algebra similar to that of Smith and Graser (1986), it can be shown that quantities arising in eqns (11) and (12) can be written in terms of (A* + >'1)-1 and For example: A! can be found using a sequence of Householder transformations, using complex arithmetic. In the Appendix, it is shown how the special structure of A l with real quadrants (L'BL and D) and imaginary quadrants (iL'C) and (iL'C)' can be used to give an algorithm for A* and q * and to evaluate eqns(12) to (16), using only real arithmetic. It should be noted that this computational strategy is using an existing iterative scheme and is manipulating the equations to reduce the number of computations.

NUMERICAL EXAMPLE
To illustrate the formation of the model and mixed model equations, we use an example of 5 observations of the same sex with individuals 2 and 3, the offspring of individual 1, and individuals 4 and 5, the offspring of individual 2. There is 1 fixed effect with 2 levels, the first 3 observations at 1 level, and the last 2 at the second level. The model (4) is then of the form: with the observation vector coded so that observations on parents (1 and 2) occur first. The variance of e i is a 5 + 3/4<r!, as animals have only one known parent. The variance of (y¡) = QA -3/4<7! + a5 + 3/4 Q A = a5 + !A, var( y s) _ 1/4 QA + QE + 3/4 QA = QE + QA . and eqn(5) becomes: Eliminating 6 1 and b 2 gives: This shows the structure of eqn(6) with two real quadrants and two imaginary quadrants in the matrix on the left hand side. For this example, estimates of the effects can be found by partitioning the coefficient matrix into 2 x 2 real and imaginary parts and using results on inverses of partitioned matrices.

DISCUSSION
This is a novel approach with the advantage of working with matrices of size 2p, rather than (n+p), and the computations involved, are of the order of (2p) 3 , rather than (n + p) 3 . This technique is, therefore, of more use in populations with a low proportion of animals used as parents. The technique can be easily extended to estimates 2 multivariate residual and additive variance components matrices when measurements are taken on all animals using the procedures developed by Meyer (1985).
Two assumptions are made. First, that all non-parents have equal parental information. If this assumption is not satisfied, non-parents with unknown parents can have dummy parents inserted into the model. In most animal breedings sets where inbreeding is consciously avoided, the second assumption (parents are not inbred) is unlikely to be important. There is more likely to be concern about residual variance homogeneity, than about inbreeding generated genetic variance homogeneity. Graser et al (1987), suggested a derivative-free method of estimating variance components, based on sequentially calculating the likelihood. Their method is an obvious competitor, but it is difficult to say precisely when each method is to be preferred. The time required, depends on the number of animals, structure of population the sequence used in calculating the likelihood, the number of variates measured and the speed of convergence of the iterative procedure. Their comments suggest that our method could be advantageous if the number of parents is less than 350.
The method has been presented for a model with additive genetic covariances, but in many cases, other components, such as litter variances, need estimation. In such cases, the procedure can be used for a given ratio of litter variance to residual variance and repeated for different values of the ratio, in a similar manner to that suggested by Smith and Graser (1986).

Reduction of complex matrices to tridiagonal form
A square matrix is said to be tridiagonal (in the first t rows) if the only non-zero elements are in the r, s elements, where s = r y-1, r, or r + 1 (r < t). A real symmetric matrix A i , of size n x n, can be reduced to tridiagonal form An-1 by a sequence of (n &mdash; 2) Householder transformations (for example, Wilkinson and Reinsch, 1971).
In this sequence of Householder transformations, at the t th stage, the transformation P t is chosen to make the tt!' row of At+, contain all zero elements, except possibly those in the (t &mdash; 1), t, and (t + 1) columns. This operation will be called pivoting on the t th row. The non-zero elements in the j th (1 < j < n) row of A n -i are in the (j -1)t' and (j + 1) t &dquo; columns, indicating the previous (j -1) and next ( j + 1 ) pivots.
For the algorithm presented in this paper, a complex matrix A l of the form: where B, C, D are real matrices [of size (n, x n l ), (n, x n 2 ) and (n 2 x n 2 )] is to be tridiagonalized.
In order to illustrate the method, a numerical example will be given.
The matrix A l satisfying eqn(A1) with: will be used, with n, = n 2 = 4 (for convenience, when matrices are symmetrical, only the lower triangular part is given). The sequence of Householder transformations derived for real matrices could be used but, as a square root of a sum of squares, at, is used and as this could be negative, this would involved complex arithmetic, which is computationally considerably more demanding than real arithmetic. A modification of the tridiagonalization is now presented which avoids complex calculations.
There are two stages involved. In the first, a sequence of Householder transformation is used that gives At+, = P t A t Pt, where At+, and P t are found to have the same form as A in eqn(A1) with quadrants of real and imaginary numbers. This involves changing the ordering of the pivoting. The rows of A l are split into two sets 1 to n l and (1 + n i ) to (n, + n 2 ), corresponding to the division into quadrants in eqn(Al). The pivoting is started on row 1 and continues in the first set until at < 0, then we pivot on the first row of the second set and continue in this set until again at < 0. The process is then repeated until nl + n2 -1 = n -1 rows have been pivoted on. This procedure produces a matrix At+, with as many nonzero elements as the usual procedure, but A t+i is not necessarily tridiagonal. For example, for the numerical example, the ordering of pivoting is found to be rows 1, 2, 5, 3, 6, 7, giving A 7 .
There are atmost, 3 non-zero elements in each row of A 7 . In the second stage, A n -i will be permuted to give a tridiagonal matrix.
The first stage is now illustrated using recursive arguments. Some housekeeping notation is needed to identify pivoted rows. In the real case, the index t was related to the number of transformations carried out (t &mdash; 1), the rows which have been transformed (1... t -1) and the next row (t) to be pivoted. The complex case is more complicated; for example, formation of A 4 for the numerical example involves pivoting on row 1, 2 and 5 in turn, and the next row to be pivoted is row 3. It is convenient to define R t , S t and T t to indicate that previous operations have used the first (R t -1) rows of the first set of rows and the first (S t -1) rows of the second set of rows as pivots and that T t is the next pivot. K t is used to indicate if the next pivot (T t ) is in the first set of rows (K t = 1) or in the second set of rows (K t = 2). Within the two sets of rows, the rows are transformed in sequence so that T t = R t , if K t = 1, and T t = n 1 + S t , if K t = 2. Hence, in the numerical example, R 4 = 2 + 1; S 4 = 5 -4 + 1; T 4 = 3; K 4 = 1; R 5 = 4. As a total of (t -1) rows of A l have been pivoted, this equals the sum of rows pivoted in the two sets; i e, (R t -1) + (S t -1) = t -1. Initially, t = 1, R l = 1, S l = 1 and the first row in the first set can be used as the first pivot, and so, T i = 1, K l r 1.
Two slightly different strategies are needed for K t = 1 and K t = 2. The derivation of At+1 is given first when K t = 1. To simplify notation in this t th stage, let r = Rt, s = S t , k = K t+1 , and suppose At is of the form (Al). As K t = 1, then the next row to be pivoted is r, so T t = r; the number of rows in the two sets already used as pivots in At+, will be r, and s -1, so that R t+1 = r + 1 and S t + 1 = s.