WOMBAT—A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML)

WOMBAT—A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML) first-author †E-mail:kmeyer@didgeridoo.une.edu.au †E-mail:kmeyer@didgeridoo.une.edu.au kmeyer@didgeridoo.une.edu.au Animal Genetics and Breeding Unit, University of New England, Armidale, New South Wales 2351, Australia Animal Genetics and Breeding Unit, University of New England, Armidale, New South Wales 2351, Australia WOMBAT is a software package for quantitative genetic analyses of continuous traits, fitting a linear, mixed model; estimates of covariance components and the resulting genetic parameters are obtained by restricted maximum likelihood. A wide range of models, comprising numerous traits, multiple fixed and random effects, selected genetic covariance structures, random regression models and reduced rank estimation are accommodated. WOMBAT employs up-to-date numerical and computational methods. Together with the use of efficient compilers, this generates fast executable programs, suitable for large scale analyses. Use of WOMBAT is illustrated for a bivariate analysis. The package consists of the executable program, available for LINUX and WINDOWS environments, manual and a set of worked example, and can be downloaded free of charge from http://agbu.une.edu.au/~kmeyer/wombat.html http://agbu.une.edu.au/~kmeyer/wombat.html Software Variance components Genetic parameters Mixed model Restricted maximum likelihood INTRODUCTION 1 One of the fundamental tasks in quantitative genetics involves the estimation of (co)variance components, i.e., partitioning of the total, observed variation into its causal components, in particular variation due to genetic effects. This provides the basis for classical estimation of genetic parameters such as heritabilities and genetic correlations as well as, more recently, applications such as mapping of quantitative trait loci (QTL). Today, variance component analyses for continuous traits are predominantly based on the mixed model, and use maximum likelihood or related methods of inference. In particular, the so-called residual or restricted maximum likelihood (REML) is widely used for analyses of continuous traits. Impetus for this was provided by applications from livestock improvement. Hence the ‘standard’ linear model, relating observations for an individual to its genotype, is commonly referred to as the ‘animal model’. One of the main attractions of the ‘animal model’ for quantitative analyses is that it allows information from all types of relatives in complex pedigrees to be utilized simultaneously. Uptake of mixed model analyses by REML has been furthered by the availability of appropriate software. While implemented in general statistical packages such as SAS, Genstat or R, more specialized packages such as ASREML (Gilmour et al., 2006 2006 ), VCE (Kovacs and Groeneveld, 2003 2003 ) and DMU (Madsen and Jensen, 2006 2006 ), tailored to quantitative models are widely used. One of the earliest of such packages available was DFREML (Meyer, 1988 1988 ). This paper describes its successor, named WOMBAT. FEATURES 1 WOMBAT is a single program, written in FORTRAN95. Its main purpose is the estimation of covariance components and the resulting genetic parameters for continuous traits, fitting a linear, mixed model. WOMBAT has been developed with the analysis of data from animal breeding programs in mind, and accommodates most models commonly fitted for such data. However, it is readily applicable to similar problems in other areas of applied statistics. WOMBAT allows for a wide range of models, accommodating numerous traits, multiple fixed and random effects, selected genetic covariance structures, random regression models and reduced rank estimation. WOMBAT employs up-to-date methods of computation, teamed with the use of efficient, commercial compilers to generate ‘fast’ executable programs. It is well suited to analyses of large data sets fitting relatively simple models. Models and types of analyses 2 WOMBAT assumes a mixed, linear model of analysis. Fixed effects accommodated range from simple means to multiple, cross-classified or nested effects and covariables. Similarly, multiple random effects with a choice of covariance structure can be fitted. By default, WOMBAT treats random effects as distributed proportionally to an identity matrix or, if flagged appropriately, proportionally to the numerator relationship matrix among individuals. If applicable, inbreeding coefficients are calculated and the inverse of this matrix is set up from a list of pedigree information, for either the animal or a sire model. Other general structures are readily fitted: WOMBAT allows for multiple, user-defined covariance structures for which the corresponding inverse can be supplied. Examples where this is appropriate include the estimation of sex-linked genetic variances or analyses fitting marker or QTL effects assumed to be distributed proportionally to their ‘Identity by Descent’ matrix. WOMBAT accommodates standard uni- and multivariate analyses as well as analyses of ‘functional’ data, fitting a so-called random regression (RR) model. For the latter, a choice of basis functions is available: Legendre polynomials and uniform B-splines are provided as in-built options. Alternatively, any function chosen and specified by the user can be used. By default, covariance matrices are estimated as unstructured: while constrained to be positive semi-definite, they are assumed to have full rank and, for q traits or RR coefficients, represented by q ( q +1)/2 parameters. Recently, there has been increasing interest in reduced rank estimation and factor-analytic models to model covariance matrices more parsimoniously. WOMBAT provides the option of modeling a covariance matrix through its leading r < q principal components only, reducing the number of parameters to r (2 q − r +1)/2. WOMBAT implements a number of special features: Different models of analysis for different traits are readily fitted, and analyses combining traits recorded repeatedly with traits with a single record only are feasible. Random effects with the same covariance structure, such as direct and maternal additive genetic effects, can be correlated. Options are provided to let the user determine which levels of fixed effects to be set to zero to account for rank deficiencies in the incidence matrix. This can reduce numerical problems, in particular for analyses with nested effects. Similarly, levels of random, genetic effects can be selected to be treated as fixed. This helps, for instance, in reducing bias when base animals have been subject to selection or allows for genetic groups to be fitted. Furthermore, the user can define additional functions of the covariance components estimated to be evaluated, together with their approximate sampling errors. Typical examples include the repeatability or the total heritability in analyses fitting maternal genetic effects. Other tasks performed 2 In addition to estimation of covariance components, WOMBAT can perform several auxiliary tasks: 1. Simple best linear unbiased prediction (BLUP) for random effects, and generalized least squares estimation of fixed effects fitted. 2. Observations for a given data set can be replaced with simulated data, sampling the random effects specified in the model of analysis (parameter file) from multivariate normal distributions. 3. Results from analyses of subsets of traits can be pooled into positive semi-definite matrices, using the method of ‘iterative summation of expanded part matrices’ (Mäntysaari, 1999 1999 ). Optionally, this can be carried out weighing individual analyses differentially. User interface 2 The user interface for WOMBAT consists of a parameter file and a number of run-time, command line options. The parameter file supplies all information on the data and pedigree files, and the model of analysis to be fitted. In addition, this file specifies the initial guess (‘starting values’) for the covariance components to be estimated and, where applicable, extra constraints to be applied or genetic parameters to be calculated. It needs to be set up following a set of simple but specific rules. These are described in detail in the ‘User Notes’ accompanying the program. Command line options are available to modify the default behavior of WOMBAT for specific analyses. These include basic options to specify a continuation run, regulate the level of screen output or select an auxiliary task to be carried out. In addition, command line options allow the user to choose the strategy used to order equations in the mixed model, to select the algorithm used to search for the maximum of the likelihood function, to modify the convergence criterion or maximum number of iterates allowed, or to enforce a specific parameterization. Technical details 2 WOMBAT incorporates state-of-the-art numerical procedures to ensure both computational efficiency and reliable maximization of the likelihood function. To reduce memory requirements and thus facilitate analyses of large data sets and models, WOMBAT employs sparse matrix storage and manipulation for the mixed model equations. In addition, most arrays are allocated dynamically to the exact size required for a particular analysis. This minimizes the memory used and, in essence, imposes limits on size of the mixed model which can be analyzed only through the amount of memory (RAM) which can be accessed. The order in which equations in the mixed model are processed can have a dramatic impact on the computational requirements, both memory and computing time required, per likelihood evaluation or REML iterate. WOMBAT employs several specialised, public domain routines to achieve orderings of the mixed model equations, which minimize the ‘fill-in’ generated and thus the operation counts of sparse matrix factorization or inversion. In particular, METIS (Karypis and Kumar, 1998 1998 ) has been found to be well capable of generating high quality orderings for large scale analyses of data from livestock improvement schemes (Meyer, 2005a 2005a ). WOMBAT utilizes the fast algorithm of Tier ( 1990 1990 ) to calculate inbreeding coefficients. It is thus capable of setting up the numerator relationship matrix between individuals for very large pedigrees quickly. The default algorithm used by WOMBAT to locate the maximum of the likelihood function is the so-called ‘average information’ algorithm, developed by Thompson and co-workers (Thompson et al., 2005 2005 ). This is implemented using an automatic differentiation (Smith, 1995 1995 ) rather than sparse matrix inversion of the coefficient matrix in the mixed model equations. To ensure an increase in the likelihood in each iterate, step sizes are scaled if necessary, using the backtracking line search of Dennis and Schnabel ( 1996 1996 ) to determine the optimum scale factor. In addition, the average information matrix is forced to be safely positive definite by adding a small term to its diagonal if a small eigenvalue or large condition number is encountered. Together, these measures ensure reliable converge for a wide range of analyses, even for difficult problems involving numerous traits or multiple random effects. In addition, WOMBAT implements both the standard and ‘parameter expanded’ (Liu et al., 1998 1998 ) expectation maximization algorithm. These algorithms are most useful in the first few iterates of an analysis, in particular when adequate starting values for the covariance components are not available, and scenarios where the average information algorithm fails. They require the sparse matrix inverse (Misztal and Perez-Encisco, 1993 1993 ) of the coefficient matrix in each iterate and, often, a large number of iterates. Furthermore, WOMBAT incorporates two derivative-free search strategies, ported from DFREML. Approximate, lower bound sampling errors of parameters and covariance components are derived at convergence from the inverse of the average information matrix. For quantities which are not linear functions of the parameters estimated, e.g. heritabilities or correlations, sampling errors are derived by first approximating the function by its first order Taylor series expansion, and then computing the variance of the latter. Applications 2 WOMBAT has been released in mid-2006. Hence there are few published analyses using this software so far. WOMBAT has been used to estimate genetic parameters for large data sets on beef cattle. This involved up 14 traits in multi-trait analyses, considering a number of reduced rank and factor-analytic models (Meyer, 2005b 2005b ; 2007 2007 ). EXAMPLE 1 In the following, the input required and output provided by WOMBAT are illustrated for a simple, bivariate analysis of data from a selection experiment in mice (Meyer, 1991 1991 ), fitting a model with three cross-classified fixed effects and additive genetic effects and permanent environmental effects of the litter as random effects. The data and files shown are available as Example 2 in the suite of worked examples for WOMBAT. Parameter file 2 The parameter file for the example is shown in Fig. 1 1 . Specific ‘codes’—shown in capital letters—are employed to identify the kind of information provided. Some information is given on the same line as the ‘code’. Other, more complex information is given in a block comprising several lines. In this case, ‘codes’ are used to mark the beginning and end of a block. As far as possible, ‘codes’ are simple abbreviations of corresponding English key words. For instance: 1. COM specifies a COMMENT on the analysis, to be printed on the output files. 2. ANAL gives the type of ANALYSIS to be carried out, which is a MULTIVARIATE analysis of 2 traits (MUV 2) in this case. 3. PEDS gives the name of the PEDIGREE file. 4. DATA gives the name of the DATA file and indicates the start of a block of lines describing the layout of the file. For each of the two TRAITS (TR1 or TR2), this block contains one line for each column. The line gives the name by which this column is referred to subsequently and, if the column represents a fixed or random effect fitted, the maximum number of its levels. For additive genetic effects (animal in this case), this is determined from the list of pedigree information. Hence a ‘0’ has been given. The block is closed by a line beginning with an END. 5. MODEL indicates the start of a block which specifies the MODEL of analysis. This block contains one line for each effect fitted and one line for each trait: (1) FIX for FIXED effects fitted, (2) RAN for RANDOM effects fitted, and (3) TR for a TRAIT to be analysed, followed by the trait number. A genetic effect, assumed to have covariance matrix proportional to the NUMERATOR RELATIONSHIP MATRIX is flagged by the qualifier NRM. Again the block is terminated by an END statement. 6. VAR followed by the name of a random effect (or ‘residuals’) and the size of the covariance matrix, indicates that the following line(s) give the starting values for the respective (co)VARIANCE components. Only the elements of the upper triangle of each covariance matrix are given. Summary of pedigree and data structure 2 WOMBAT provides files with basic summary information on the pedigree and data structure encountered. These are written after the initial ‘set-up’ steps. Performing the set-up steps only and inspecting these files can help to ensure that the model of analysis has been specified and interpreted correctly. Fig. 2 2 shows part of the pedigree summary file for the example. Given are the numbers of animals in the data and without records, the numbers of sires and dams, and parents and grand-parents (not shown) identified. In addition, numbers of inbred animals, mean inbreeding coefficient and characteristics of the inverse of the numerator relationship matrix (or its user-supplied equivalent) are reported (not shown). Summary statistics reported for the data structure (not shown) comprise means, standard deviations and ranges for traits and covariables fitted, and numbers of records for each pair of traits. Random effects in the model are listed together with the respective numbers of levels found and the covariance structure specified. Correspondingly, fixed effects fitted and the number of levels identified in the data are listed separately for each trait. Estimates of covariance components 2 Part of the main results file for the example is shown in Fig. 3 3 . After some information identifying the particular analysis, such as the comment line from the parameter file and the input files used, the maximums of the REML log likelihood together with the corresponding Akaike and Bayesian information criteria (AIC and BIC) are given. Not shown here, this is followed by a list of the parameters as estimated, and details of the convergence criteria for the last 3 iterates. Subsequently, estimates of the covariance matrices are reported, beginning with the residual covariances and the other covariance matrices in the order in which the corresponding starting values have been specified in the parameter file. In addition, WOMBAT attempts to calculate a matrix of total, phenotypic variances as the sum of all covariance matrices estimated. For each matrix, its eigenvalues and the proportion of variance explained, and the corresponding correlations and variance ratios (i.e. proportions of the total variance) are given. In addition, the covariance components and correlations or variance ratios are listed together with their approximate lower bound sampling errors. Estimates for fixed and predictions for random effects fitted 2 WOMBAT obtains solutions, i.e. estimates for all fixed effects and predictions for all random effects fitted as a by-product, and these are written out to separate files at the end of the analysis. Fig. 4 4 shows the estimates of fixed effects for the example. For each effect level, listed are the effect name, original level code, running number, number of records for the level, and the ‘raw’ mean. Solutions are given as estimated, setting the first level of fixed effects other than the first (i.e. the second effect onwards) to zero to account for rank deficiencies, and transformed so that they sum to zero. As WOMBAT uses an implementation of the average information REML algorithm which does not involve inverting the coefficient matrix of the mixed model equations, standard errors of estimates are not provided automatically. However, options are available to obtain these if required. AVAILABILITY 1 WOMBAT is available to the scientific community free of charge. It can be obtained by downloading from the WOMBAT home page ( http://agbu.une.edu.au/~kmeyer/wombat.html http://agbu.une.edu.au/~kmeyer/wombat.html ) Material available comprises the program, documentation and a set of examples. Program 2 WOMBAT is provided as a pre-compiled (‘executable’) program. WOMBAT has been developed in a LINUX environment and is intended to be used under such operating system. It has been tested under Fedora Core 2 to 5 and Ubuntu 7.04. LINUX versions available are: 1. A 64-bit version and 32-bit version compiled using the Pathscale compiler. 2. 32-bit versions, ‘generic’ or optimized for a Pentium IV processor, compiled using the Lahey LF95 compiler. These are limited to analyses using no more than 2 Gb of RAM. 3. A 32-bit version, compiled with all debugging and checking options switched on. This is too slow for general use and intended for use with bug reports only. In addition, WINDOWS versions of WOMBAT, compiled using the g95 compiler, are available. These are intended for demonstration and evaluation of the program or relatively small scale analyses, and have thus been restricted to models with less than 150 000 effects (equations) in total. WINDOWS versions available are: 1. A 32-bit version to run under CYGWIN, 2. A 32-bit version, compiled under MinGW, to run in an MSYS or MsDOS window. Documentation 2 A manual (‘User Notes’) for WOMBAT is available as a .pdf file. Currently, the manual comprises just over 60 pages (A4). It describes the tasks performed, options available, input required and output provided. In addition, an appendix gives some technical details on the methods implemented. It is intended to provide additional comments and documentation related to WOMBAT and more general issues of REML estimation in form of a Wiki, i.e., a set of web pages. This includes information on updates or changes to the program, hints on special types of analyses, and general comments or links. While the main Wiki pages are restricted, there is a freely editable ‘open’ section available. The Wiki can be accessed through a link from the WOMBAT home page. It represents work in progress; to date contents are limited. Examples 2 A suite of worked examples is provided. Examples comprise data, pedigree and parameter files (input), as well as output files generated by WOMBAT. Currently there are 9 examples, most with multiple sub-analyses. By and large examples represent practical types of analyses common to animal breeding applications. Features illustrated include multivariate analyses, random regression models, analyses fitting multiple random effects, different models for different traits, and reduced rank estimation. References JE Dennis , RB Schnabel . Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SIAM Classics in Applied Mathematics , 1996 . : Philadelphia Society for Industrial and Applied Mathematics , . Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SIAM Classics in Applied Mathematics AR Gilmour , DJ Gogel , BR Cullis , 等 . ASReml User Guide Release 2.0 , 2006 . : Hemel Hempstead, HP1 1ES, U VSN International Ltd. , 342 . ASReml User Guide Release 2.0 G Karypis , V Kumar . METIS, a Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-in Reducing Ordering of Sparse Matrices Version 4.0 , 1998 . : Minneapolis Department of Computer Science, University of Minnesota , 44 . METIS, a Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-in Reducing Ordering of Sparse Matrices Version 4.0 MA Kovacs , E Groeneveld . VCE-5 User’s Guide and Reference Manual, Version 5.1 , 2003 . : Mariensee,Germany Institute of Animal Husbandry and Animal Behaviour, Federal Research Center of Agriculture , 68 . VCE-5 User’s Guide and Reference Manual, Version 5.1 C Liu , DB Rubin , YN Wu . Parameter expansions to accelerate EM: the PX-EM algorithm . Biometrika , 1998 . 85 ( 4 ): 755 - 770 . DOI: 10.1093/biomet/85.4.755 http://doi.org/10.1093/biomet/85.4.755 . Parameter expansions to accelerate EM: the PX-EM algorithm P Madsen , J Jensen . A User’s Guide to DMU: A Package for Analyzing Multivariate Mixed Models, Version 6 , 2006 . : 8830 Tjele, Denmark Danish Institute of Agricultural Sciences, Dept. of Animal Breeding and Genetics, Research Centre Foulum , 27 . A User’s Guide to DMU: A Package for Analyzing Multivariate Mixed Models, Version 6 EA Mäntysaari . Derivation of Multiple Trait Reduced Random Regression (RR) Model for the First Lactation Test Day Records of Milk, Protein and Fat . 50th Annual Meeting , 1999 . Europ. Ass. Anim. Prod. Mimeo , 8 . Derivation of Multiple Trait Reduced Random Regression (RR) Model for the First Lactation Test Day Records of Milk, Protein and Fat K Meyer . DFREML—a set of programs to estimate variance components under an individual animal model . J Dairy Sci , 1988 . 71 ( Suppl. 2 ): 33 - 34 . . DFREML—a set of programs to estimate variance components under an individual animal model K Meyer . Estimating variances and covariances for multivariate animal models by restricted maximum likelihood . Genet Sel Evol , 1991 . 23 ( 1 ): 67 - 83 . DOI: 10.1051/gse:19910106 http://doi.org/10.1051/gse:19910106 . Estimating variances and covariances for multivariate animal models by restricted maximum likelihood K Meyer . Ordering strategies to reduce computational requirements in variance component estimation . Proc Ass Advan Anim Breed Genet , 2005 . 16 282 - 285 . . Ordering strategies to reduce computational requirements in variance component estimation K Meyer . Genetic principal components for live ultra-sound scan traits of Angus cattle . Anim Sci , 2005 . 81 ( 3 ): 337 - 345 . DOI: 10.1079/ASC50850337 http://doi.org/10.1079/ASC50850337 . Genetic principal components for live ultra-sound scan traits of Angus cattle K Meyer . Multivariate analyses of carcass traits for Angus cattle fitting reduced rank and factor-analytic models . J Anim Breed Genet , 2007 . 124 ( 2 ): 50 - 64 . DOI: 10.1111/j.1439-0388.2007.00637.x http://doi.org/10.1111/j.1439-0388.2007.00637.x . Multivariate analyses of carcass traits for Angus cattle fitting reduced rank and factor-analytic models I Misztal , M Perez-Encisco . Sparse matrix inversion for restricted maximum likelihood estimation of variance components by Expectation-Maximization . J Dairy Sci , 1993 . 76 ( 5 ): 1479 - 1483 . . Sparse matrix inversion for restricted maximum likelihood estimation of variance components by Expectation-Maximization SP Smith . Differentiation of the Cholesky algorithm . J Comp Graph Stat , 1995 . 4 ( 2 ): 134 - 147 . DOI: 10.2307/1390762 http://doi.org/10.2307/1390762 . Differentiation of the Cholesky algorithm R Thompson , S Brotherstone , IMS White . Estimation of quantitative genetic parameters . Phil Trans Roy Soc B , 2005 . 360 ( 1459 ): 1469 - 1477 . DOI: 10.1098/rstb.2005.1676 http://doi.org/10.1098/rstb.2005.1676 . Estimation of quantitative genetic parameters B Tier . Computing inbreeding coefficients quickly . Genet Sel Evol , 1990 . 22 ( 4 ): 419 - 425 . DOI: 10.1051/gse:19900404 http://doi.org/10.1051/gse:19900404 . Computing inbreeding coefficients quickly 10.1631/jzus.2007.B0815 Project (No. BFGEN.100B) supported by the Meat and Livestock Ltd., Australia (MLA) 2007-09-12 2007-09-14 2022-01-05 Copyright © 2007, Journal of Zhejiang University Science 2007 Journal of Zhejiang University-SCIENCE B (Biomedicine and Biotechnology) 11 08