Jaffrezic, F.

Thompson, R.

Rothamsted Repository

Patron:		Her	Majesty	The	Queen	 	 Rothamsted	Research	Harpenden,	Herts,	AL5	2JQ		Telephone:	+44	(0)1582	763133	Web:	http://www.rothamsted.ac.uk/			 			Rothamsted Research is a Company Limited by Guarantee Registered Office: as above.  Registered in England No. 2393175. Registered Charity No. 802038.  VAT No. 197 4201 51. Founded in 1843 by John Bennet Lawes.		Rothamsted Repository DownloadC1 - Edited contributions to conferences/learned societiesThompson, R. and Jaffrezic, F. 2003. Modelling and estimation for the genetic analysis of longitudinal data. The publisher's version can be accessed at:• http://isi.cbs.nl/iamamember/CD3/index.htmlThe output can be accessed at: https://repository.rothamsted.ac.uk/item/8922v.© 14 August 2003, Please contact library@rothamsted.ac.uk for copyright queries.17/07/2019 11:40 repository.rothamsted.ac.uk library@rothamsted.ac.ukModelling and Estimation for the Genetic Analysis ofLongitudinal DataRobin Thompson 1,2 and Florence Jaffrezic  31 Rothamsted Research, Harpenden,Herts AL5 2JQ,England, 2Roslin Institute (Edinburgh), Roslin,Midlothian EH25 9PS,Scotland,3 INRA,SGQA,78352,Jouy-en-Josas,France.1     INTRODUCTIONIn animal breeding situations we are interested in predicting future performance .Over time thenumber and types of data available have gradually increased .We have expanded from analyseson  single traits  ,taking into account relationships of animals with their father and ignoringenvironmental effects into analyses that include several traits ,taking account of all geneticrelationships and realistic  environmental models .The basic model  is a mixed linear modelwith predictions based on best linear unbiased prediction (BLUP,Henderson,1973).Residualmaximum likelihood (REML,Patterson and Thompson ,1971) is often the method of choice forestimation of variance parameters.We discuss algorithms that allow estimation in lage indusrialsized populations.There has recently interest in developing estimation methods for longitudinaldata ,for example the so called test day models for milk yield in dairy cattle where informationof a series of measurements are combined to give predictions for total lactation yield and othercomponents of the lactation curve.2  RESIDUAL MAXIMUM LIKELIHOOD ESTIMATION.We consider a linear model  eZuXby ++=  with var(y)= ZGZ‘ +R , var(u)= G and var(e)=R.The residual log-likelihood(REML) is of the formX)1VXlogdet(logdet(V))bˆX(y1V)bˆX(yαL −′−−−−′−This is different from the usual likelihood form in that it is a function of error contrasts –contrasts that do not tell us about fixed effects.  This difference has two consequences, the useof the weighted least squares estimate of b, bˆ ,given by  y1VXbˆX1VX −′=−′The term in det( X1VX −′ )that is sometimes thought of as a penalty function because the fixedeffects are not known.  Mixed model equations (Henderson, 1973) pay an important part in theanalysis process.  These are of the form −′−′=−+−−−−y1RZy1RXuˆbˆ1GZ1RZ'X1RZ'Z1RX'X1RX'Terms derived from these include prediction error variances found from writing the mixedmodel equations as   Cs = R  so that1Cuuˆbˆvar −=−It is often useful to express relevant quantities in terms of the projection matrix1VX1X)1VXX(1VP −′−−′−−=X)1VXlog(logdet(V)PyyαL −′−−′Estimation of a variance parameter θi involves setting to zero the first derivatives[ ])iθV/P(tr)PyiθV/P(yiθL/ ∂∂−∂∂′=∂∂These could be thought of equating a function of the data to its expectation. normally finding amaximum of the likelihood requires an iterative scheme. One suggested by Patterson andThompson (1971) this is based on the expected value of the second differential that is)P]iθV/)P(iθV/tr[P((1/2))2iθL/2E( ∂∂∂∂−=∂∂This is called the Expected Information.  Using the first and second differentials we can updateθ using the rate that all the terms from solution of MME and C-1 for exampleθ).L/(1EInfθθˆ ∂∂−+=Whilst this development is very direct, later developments tried to take account of the structureto reduce the computational effort An alternative algorithm was suggested by Dempster, Lairdand Rubin (1987).This EM algorithm is based on thinking of the random effects as `missing’.The estimation is based on using  PEV(u)uu2gσˆs +′=  writing this as],)iθG/(1tr[V2gsσy1)ViθG/(1Vy2gσˆs ∂∂−−+−∂∂−′=We see this as a manipulation of equating the first differential to zero.    It can be also writtenθ)L/(1Infθθˆ ∂∂−+=  with Inf representing the information on the complete data.  Oneadvantage of this method is  2gσ  that stays in the parameter space 02gσ ≥ .Another advantage is that there is an increase in likelihood in each iteration.  Disadvantages arethat the method can be slow to converge (indeed this method is said to be the most widely usedin terms of numbers of iterations) and it requires the inversion of C in each iteration.An important advance was the rediscovery (Misztal and Perez-Enriso, 1993) of an algorithm(Takahashi, et al. 1973) that allowed the calculation of the `relevant’ terms in the inverse of Crequired for forming the first differentials without calculating all the elements of the inverse.This result allowed the implementation of EM algorithms to estimate variance parameters,(Misztal, 1994) for bigger problems.  These were an improvement on derivative free methodsbut could still be slow to converge. To reduce the computation of the information matricesThompson and co-workers (Johnson and Thompson, 1995, Gilmour et al., 1995, and Jensen etal., 1997) suggested using an alternative information matrix. The second differential of C with respect to θi and θj.)]jθV/)P(iθV/(1/2)tr[P()]jθiθL/2E[(andPy)jθV/)P(iθV/P(y)]jθV/)P(iθV/(1/2)tr[P()jθiθL/2(∂∂∂∂−=∂∂∂∂∂∂∂′−∂∂∂∂=∂∂∂ Both these terms often called observed and expected information are difficult to calculate butthe average     Py)jθV/)P(iθV/P(y(1/2))]jθiθL/2AI[([ ∂∂∂∂′−=∂∂∂    can be calculated byusing )PyiθV/( ∂∂ and Py)jθV/( ∂∂ as working variables and obtaining the residual cross-product between these working variables.  This calculation is much simpler than calculatingeither the observed and expected information. A synthesis of comparisons of these iterativemethods was carried out by Hofer (1998) .  These show the expected improvement of EMmethods over derivative free methods.  They also show that most second differential methodsconverge in relatively small number of iterations.Computer software is available  to fit verygeneral mixed models using the AI algorithm (Gilmour et al. 2002)3     LONGITUDINAL DATAWe assume the  observed pheotypic tragectory can be decomposed as Y(t)=µ(t) + g(t) + ,e(t) +ε,where µ(t) is a function of t,g(t) is the genotypic mean function of  Y(t), ε is the residualvariation ,assumed normally distributed with unknown variance,g(t) and e(t) are Gaussianvariables and represent time (t) dependent genetic and enviromental deviations ,independentod one another with covariance functions G(s,t) and E(s,t). Basically three types of modelshave been suggested for dealing with genetic longitudinal data.3.1 Random Regression models. Random regression (RR)models are well known in the context oflongitudinal data analysis(Diggle et al.,1994).Often convenient  parametric curves as linearfunctions of t are chosen to summerise  the genetic and enviromental deviations.By allowingvariances and covariances amongst the regression coefficients one can generate variancematrices as a function of  t.Alternatively these matrices can be generated as a way of smoothingpreviously estimated covariance matrices (Kirkpatrick and Heckman,1989,Meyer andHill,1997)3.2 Structured antedependence models. In contrast to starting with a linear formulation one mightstart with a time series like formulation and think that observations at time t might be explainedin terms of previous ones..An antedependence structure of order r is defined by the fact that theith observation (i> r) given the previous r  preceding ones is independent of all previous ones(Gabriel,1962).One can easily generalize this to deal with genetic and enviromentalcomponents.In structured  antedependance (SAD) models the conditional variuance at time t ismodelled using a parametric function of t ,for example a polynomial in t (Nunez-Anton andZimmerman,2000).3.3  Character Process models. In these models rather than thinking of a linear model  to generatethe covariance functions we directly decompose the variance function G(s,t) into terms such asυ(s)υ(t) ρ(‌s-t‌).The variance function υ(t) υ(t) can be written as a polynomial in t and ρ(‌s-t ‌)describes how the correlation changes with s and t. Jaffrezic and Pletcher(2000) suggest using atransition of t,based on a small number of parameters, as a way of introducing non-stationarity.4     COMPARISONS FOR UNIVARIATE ANALYSISF. Jaffrezic (Ph. D thesis,2001) carried out an extensive investigation of a variety of simulatedcovariance structures and empirical data and found under most circumstances  that structuredantedependence and character processes provide the best description of the underlyingcovariance structure.One example used by Jaffrezic et al.(2002)  was based on an analysis of  adata set of 9277 cows from 464 bulls with 10 records  per animal.Table1 shows that a SADmdel of first order for the genetic part and third order for the enviromentasl part has a higherorder likelihood than a quartic random regression with far fewer parameters (11 instead of 31)5 DISCUSSIONThese analyses show that CP and SAD models offer advantages to random regressions to fitcovariance structures that occur in genetic analysis with fewer parameters.We have recentlyfound that such advantages carry over to multivariate analyses(Jaffrezic et al., ,2003) buildingon work of Sy et al.(1997) with even bigger reductions of parameters from RR models andoften very similar models being useful for different traits. One advantage of RR models is thatfor half-sib data sire and animal models  are equivalent.This simplification does not hold for CPmodels.In one particular case Jaffrezic et al. (2002) the animal model parameterization wasbetter  because it allowed a more appropriate model for the enviromental  covariance matrix.6     CONCLUSION        We have discussed estimation methods and models that allow more flexible models of geneticdata.REFERENCESDiggle, PJ, Liang, KY & Zeger, SL (1994) “Analysis of Longitudinal Data”. Oxford ,ClarendonPress.Dempster, A.P., Laird, N.M. & Rubin, D.B. (1977) J. Royal Stat. Soc., Series B 39: 1-38.Gabriel,K.R.(1962) Ann. Math. Stat. 33:201-212.Jaffrezic, F. and Pletcher, S.D. (2000) Genetics 156 ,913-922.Jaffrezic, F,White IMS, Thompson R, Visscher PM (2002) J. Dairy Sci.,85,968-975.Jaffrezic. F, Thompson R. ,Hill, W.G. (2003) Genet. Res. (in press) Gilmour, A.R., Thompson, R. & Cullis, B.R. (1995)  Biometrics 51: 1440-1450.Gilmour, A.R.,Gogel, B.J.,Cullis, B.R.,Welham, S.J.,Thompson, R. (2002) ASREML User Guide1.0 .VSN International  Ltd. ,Hemel Hempsted HP1 1ES.Henderson, C.R. (1973) In:  Proc. Animal Breeding and Genetics Symposium in Honour of Dr. JayL. Lush,  Champaign, Illinois:10-41.Henderson,C.R, (1975) Biometrics 31: 3-447.Hofer, A. (1998) J. Anim. Breed. Genet. 115: 247-265.Jensen, J , Mantysaari, E , Madsen, P.,  & Thompson,R. (1997) J. of Indian Soc. Of Agric.Sci. 49: 215-236.Johnson, D.L. & Thompson, R. (1995) J. Dairy Sci , 78,: 449-456.Kirkpatrick,,m. and Heckman, N.(1989) J. Math. Biol. 27,427-450.Meyer, K. and Hill, W.G.(1997) Livest. Prod. Sci. 149, 185-200.Misztal, I. & Perez-Enciso, M. (1993) J. Dairy Sci. 76: 1479-1483.Misztal, I. (1994) J. Anim. Breed. Genet. 111: 346-355.Nunez-Anton, V. and Zimmerman, D. I.(2000) Biometrics, 56,699-705.Patterson, H.D. & Thompson, R. (1971). Biometrika, 58:545-554Smith, S.P. & Graser, H.-U. (1986) J. Dairy Sci. 69: 1156-1165.Sy ,J.P.,Taylor,J.M.G.,Cumberland,W.G.(1997) Biometrics , 53:542-555.Takahashi, K., Fagan, J. & Chin, M.S. (1973) In  Proc. 8th Inst. PICA Conf.Minneapolis: 63.RÉSUMÉLe but de cet exposé est tout d’abord de présenter les différentes procédures d’estimation proposéespour l’estimation REML des composantes de la variance, dont l’algorithme EM et un algorithme desecond ordre basé sur la matrice d’Average Information. Nous présenterons ensuite trois types demodèles proposés pour l’analyse des données longitudinales : la aléatoire, les modèlesantédépendants structuraux et les modèles à processus. La comparaison de ces différents approchesdans une analyse génétique de données de production laitière montre l’intérêt des modèlesantédépendants qui permettent une meilleure modélisation de la structure de covariance que lesmodèles de régression aléatoire, tout en nécessitant beaucoup moins de paramètres.Table 1 Model comparisons for the genetic analysis of lactation curves for dairy cattle(NPCov:number of parameters in the covariance structure , LogL :Log Likelihood,US Unstructuredcovariance matrix,SAD(I): structured model of order I)Model Genetic Environment NPCov LogLUnstructured1 US US 110 41262 SAD(1) US  59 4109Structured Antedependent3 SAD(1) SAD(3)  11 38454 SAD(2) SAD(3)  12 38525 SAD(3) SAD(3)  13 38546 SAD(2) SAD(2)  11 37967 SAD(1) SAD(1)   9 3580 Random Regression8 Quartic Quartic  31 36239 Quadratic Quartic  22 360710 Cubic Cubic  21 333611 Quadratic Quadratic  13 2767

Modelling and estimation for the genetic analysis of longitudinal data

https://repository.rothamsted.ac.uk/download/04f2a65bf03827b4eea54e6317b383d431246ae603a440655da50c5e8a62e6e3/3544454/3490.pdf

Modelling and estimation for the genetic analysis of longitudinal data

Abstract

Similar works

Full text

Available Versions

Rothamsted Repository