Chapter 5
Run options

5.1 Overview

Running WOMBAT can be as simple as specifying its name at command level, i.e.

EXAMPLE:

wombat

A number of options are available, however, to modify the run time behaviour of WOMBAT. These range from simple options to set the verbosity level of screen output or to select a continuation run, to highly specialised options to modify the search strategies for the maximum of the likelihood function. Multiple options can be combined, but some care is required so that options given last do not unintentionally cancel out some of the effects of options given first. Options available are summarised in 5.1, 5.2 and 5.3.

Run options can be given in three ways :

1.

On the command line, e.g. wombat -v -c myparfile.par

2.

In a file with the default name RunOptions in the current working directory.
This file must contain a single option per line. If such file exists, it is read before any command line options are parsed. There are no checks for options specified in duplicate. This implies that the last setting for a particular option is the one to be used, and that the options in RunOptions can be overridden by command line options.

3.

On the first line of the parameter file, after a code of RUNOP (see section 4.3). Any such options are read before options specified on the command line, i.e. the latter may overide them.

Following Linux standard, run options have the form “-a” where a stands for a single, lower-case letter, or the form “--abcdef” where abcdef stands for a multi-letter code.

5.2 Basic run options

5.2.1 Continuation run

In some instances, REML estimation continuing from the ‘best’ estimates so far is necessary. This is invoked with the -c option.
If specified, WOMBAT will attempt to read its ‘starting’ values from the file BestPoint in the current working directory. N.B. If this is not available or if a read error occurs, the run proceeds as a ‘new’ run, i.e. using the starting values given in the parameter file instead. The -c option selects estimation using the AI algorithm, unless another procedure is selected explicitly.

5.2.2 Level of screen output

The amount of screen output can be regulated by the following options

-t

selects terse output

-v

selects verbose output, useful to check the parameter file

-d

selects detailed output, useful for debugging

5.2.3 Set-up steps

It is good practice, to start each analysis with a run which carries out the set-up steps only, checking that the summary information provided on the model of analysis and data structure in file SumModel.out (see subsection 7.1.2) is as expected, i.e. that WOMBAT is fitting the correct model and has read the data file correctly. This is specified using --setup. This option also invokes verbose screen output.

5.2.4 Quality of starting values

For analyses comprising 18 or less covariance components to be estimated, WOMBAT defaults to the AI algorithm for maximisation. For other analyses, the default is for WOMBAT to begin with up to 3 iterates of the PX-EM algorithm, before switching to an AI algorithm. For reduced rank estimation, every tenth AI iterate is subsequently replaced by a PX-EM step. Two options are provided to modify this according to our perception of the quality of starting values for covariance components available, without the need to consider the ‘advanced’ options below.

If we are confident, that we have good starting values, this might cause an unnecessary computational overhead. Specifying --good reduces the maximum number of (PX-)EM iterates to 1. Similarly, for potentially bad starting values, estimation might converge more reliably if a few more initial (PX-)EM iterates were carried out. Specifying --bad sets this number to 8. If the increase in log likelihood during the initial (PX-)EM iterates is less than 2.0, WOMBAT will switch to the AI algorithms immediately.

5.2.5 Numerical settings

5.2.5.1 Operational zero

WOMBAT uses a small positive value as operational zero  to reduce numerical errors. The default value is 10⁻⁸. This can be changed specifying the run option --zeron where n is a single-digit integer. This option redefines the operational zero to be 10⁻ⁿ.

5.2.5.2 Small pivots in Cholesky decomposition

For variance component estimation parameterising to the elements of the Cholesky factors of covariance matrices, estimates are constrained to the parameter space by restricting the pivots in the decomposition to a small value. The default is 0.001. This can be changed through the run option --pivotx where x is the real value specifying the new limit.

5.2.5.3 Augmenting diagonals

In some instances a run can fail because of numerical issues in the Cholesky factorisation of the coefficient matrix of the MME. Most commonly this involves LAPACK routine DPOTRF. Reasons may include a build up of rounding errors or a coefficient matrix which is not safely positive definite to begin with. A possible work around then is to add a small, positive constant to the diagonal elements of the coefficient matrix. This can be requested by adding
--fudge
to the command line (or first line in the parameter file). The default is to add 10⁻⁸. This can be changed to 10⁻ⁿ by specifying --fudgen instead with n an appropriate integer value.

5.2.6 Intermediate results

WOMBAT writes out the currently best values of estimates of covariance components to the file BestPoint whenever the likelihood is increased. The option --best causes WOMBAT to read this file and write out the matrices of covariances and corresponding correlations in more readily legible format to the file BestSoFar.out. WOMBAT with this option can be used in a directory in which an estimation run is currently active, without interfering with any of the files used.

5.2.7 Prediction only

WOMBAT can be used for a simple BLUP run, using the ‘starting’ values given in the parameter files as assumed values for the true covariances. In this mode, no pruning of pedigrees is carried out. If -c is specified in addition, WOMBAT will try to read the values from the file BestPoint instead – this is useful to obtain ‘backsolutions’ after convergence of an estimation run. Solutions can be obtained either directly or iteratively:

• Run option --blup selects the direct solution scheme. This involves sparse matrix inversion of the coefficient matrix, i.e. computational requirements are similar to those of a single REML iterate using the EM algorithm. While this can be computationally demanding for larger problems, it allows standard errors and expected accuracies (correlations between estimated and true effects) for random effects to be calculated from the diagonal elements of the direct inverse.
• Run option --solvit specifies that the mixed model equations are to be solved iteratively. The default is a preconditioned conjugate gradient (PCG) algorithm. Where applicable (multivariate or random regression models), either scheme utilises the inverse of the diagonal block comprising all equations pertaining to a particular random effects level. Up to 50000 iterates are performed. The default convergence criterion requires the square root of the sum of squared deviations in solutions between iterates divided by the sum of squared solutions to be less than 10⁻⁸. This can be overridden by specifying an integer number - to be the new exponent - immediately after the option; e.g. --solvit6 sets the convergence criterion to a less stringent value of 10⁻⁶.
  This option increases the limit on the maximum number of effects (equations) which can be fitted to a multiple of the maximum imposed for REML analyses. No attempt is made to approximate standard errors, and no checks for redundant animals in the pedigree are performed in this mode.

  HINT: For large problems, --solvit is best combined with --choozhz (see subsection 5.3.1).

• Run option --mmeout acts like --solvit, but writes out  the complete mixed model equations in a form suitable for further analyses. This includes a file containing all non-zero elements in the lower triangle of the coefficient matrix (default name MMECoeffMatrix.dat), and the file MMEEqNos+Solns.dat which specifies the mapping of equation numbers to effects fitted, as well as giving the vector of right hand sides and solutions; see subsection 7.2.9 for details.
• Run option --s1step provides an iterative solver as in --solvit, but targeted towards a so-called single-step analysis and providing a PCG algorithm only. The main difference is that equations for genotyped individuals are grouped together and are then processed as a dense block. This requires the inverse of the  combined relationship matrix (H⁻¹) to be supplied as a *.gin matrix, and a code for each individual – specifying whether is has been genotyped or not – to be given in the third, additional column of the pertaining *.codes file.

  Change: The default for the PCG algorithm has been changed to the so-called SSOR preconditioner. There are now five additional choices to vary the preconditioner by appending the letter A, B, C, D or E:

  --s1stepA

  invokes a simple diagonal preconditioning scheme.

  --s1stepB

  provides a blockdiagonal preconditioner with the blocks equal to the effects for all traits in a multivariate analysis. However, if genetic groups are fitted as an additional random effect, the diagonal block for all levels (and traits) is used.

  --s1stepC

  is like --s1stepB except that the full diagonal block for genotyped animals (all levels and traits) is utilised.

  --s1stepD

  as for --s1stepC, full blocks for genotyped animals and genetic groups are used for preconditioning, but this is combined with the diagonal preconditioner for the remaining equations

  --s1stepE

  provides a diagonal preconditioner for all effects other than genetic groups (equivalent to B for univariate nalyses).

  	A	B	C	D	E
  	’diagonals’	’diagblock’	’genoblock’	’diaggenos’	’diagggrps’
  genotyped animals	diaga	blox	full	full	diag
  non-geno. anim.s	diag	blox	blox	diag	diag
  genetic groups	diag	full	full	full	full
  other effects	diag	blox	blox	diag	diag

• Run option --s2step is similar, but uses iteration on data  to reduce the memory requirements and uses a diagonal preconditioner.
• Run option --snap selects a simple GWAS  analysis – for fixed variance components (i.e. a BLUP run)– fitting a mixed model with a single SNP effect (fitted as a linear covariable) but processing multiple SNPs with reference allele counts read in from a separate file (default name QTLAllels.dat or QTLAllelsR.dat; see subsection 6.6.3). This assumes no missing counts and employs an efficient computing strategy which exploits that only the equation for the SNP changes for different SNPs [40]. Estimated SNP effects together with their standard errors are written out to a file with the default name QTLSolutions.dat (see subsection 7.2.10). This run option must be used in conjunction with a SPECIAL statement which identifies which of the covariables in the model of analysis is to be treated as SNP effects (see subsection 4.17.7).

5.2.8 Simulation only

Specifying --simul causes WOMBAT to sample values for all random effects fitted for the data and pedigree structure given by the respective files, from a multi-variate normal distribution with a mean of zero and covariance matrix as specified in the parameter file. Again -c can be used to acquire population values from the file BestPoint instead. No fixed effects are simulated, but the overall (raw) mean for each trait found in the data set is added to the respective records. Optionally, --simul can be followed directly (i.e. no spaces) by an integer number n in the range of 1 to 999. If given, WOMBAT will generate n simulated data sets (default n = 1).

Simulation uses two integer values as seeds to initialise the pseudo-random number generator. These can be specified explicitly in a file RandomSeeds (see subsubsection 6.6.5.3). Output file(s) have the standard name(s) SimData001.dat, SimData002.dat, …, SimDatan.dat. These have the same layout as the original data file, with the trait values replaced by the simulated records; see subsection 7.2.5.

This option is not available for models involving covariance option GIN (see subsection 4.10.4).

5.2.9 Sampling to approximate standard errors

Large sample theory indicates that maximum likelihood estimates are normally distributed with covariance matrix equal to the inverse of the information matrix. Hence, sampling from this distribution has been advocated as a simple strategy to approximate standard errors of ‘complicated’ functions of variance components without the need for a linear approximation or to evaluate derivatives [35].

WOMBAT provides the run option --samplen to obtain such samples in a post-analysis step, with n denoting the number of samples to be drawn. If n is omitted, the number of samples is set to a value of 10000. It requires output files from an estimation run (at convergence), namely BestPoint, AvInfoCovs and AvInfoParms to specify the multivariate Normal distribution to sample from. By default, samples for all covariance components in the model are generated. In addition, the covariance matrix for a single random effect only can be selected by specifying its name in a SPECIAL statement (see subsection 4.17.8). Samples generated are written to a file with the standard name CovSamples_ALL.dat or CovSamples_name.dat with name the name of the random effect selected (see subsection 7.2.12). These files are ‘formatted’ (text) files that are suitable for further analysis using standard statistical software packages. In addition, a file named SumSampleAI.out is written which summarizes details of the run together with means and variances across replicates.

This option is not implemented for random regression analyses or random effects with diagonal covariance matrices. It tends to work best for estimation runs with the PC option, even when matrices are estimated at full rank.

5.2.10 Matrix inversion only

Use of this facility has been changed and the following section has been UPDATED!

WOMBAT can be used as a ‘stand-alone’ program to invert a real, symmetric matrix. Both a ‘full’ inverse and a ‘sparse’ inverse are accommodated, where the latter comprises only the elements corresponding to the non-zero entries of the Cholesky factor of the matrix given.

Input: The matrix to be inverted is expected to be positive definite (see below for exceptions) needs to be supplied in a file, with one record for non-zero element in one triangle of the symmetric matrix. Each record has to contain three space separated items:

a)

row number,

b)

column number and

c)

the matrix element.

If the matrix is to be inverted in sparse mode, it is best to supply these elements in the order of the lower triangle row-wise.

Options: There are several different ‘modes’, invoked by different run options, mostly beginning with --inv and followed (space separated) by the name of the file containing the matrix. Note that these should be at the end of the command line, i.e. that any other run options need to be given before it. No parameter file is read or expected.

1.

Using  --invert (or simply --inv) is the default option which selects inversion of a full-stored, dense matrix. This uses LAPACK routines DPOTRF and DPOTRI.

2.

Option --invlap selects the same calculations as --inv. Expanding this to --invlapp or --invlapf (note the extra letter!) changes the matrix storage to ‘packed’ form, storing only one triangle of the symmetric matrix and thus requiring less memory. For --invlapp, calculations are carried out using LAPACK routines DPPTRF and DPPTRI and for --invlapf, LAPACK routines DPFTRF and DPFTRI are used. These two differ in the in way the elements of the matrix are packed and the latter tends to be computationally faster (see the LAPACK manual for details). Note that for integer*4 addressing, the largest size of a packed matrix that can be accommodated is 65,535.

3.

Option --inveig first performs an eigen-value and -vector decomposition of the matrix, using LAPACK routine DSYEVR. If eigenvalues less than a specified value (default: 10⁻¹² ) are found, these are set to this minimum, before constructing the inverse from the eigenvectors and -values. Eigenvalues are written out to a separate file.

4.

Option --invspa chooses sparse matrix  inversion. This is suitable for the inversion of large, sparse matrices where selected elements of the inverse are sufficient. Calculated and written out are all elements of the inverse corresponding to the non-zero elements in the Cholesky factor of the original matrix; any other elements of the inverse are ignored. The methodology is the same as that used in the Expectation-Maximisation steps in REML estimation, using the supernodal implementation.

N.B.: Remember that this option does not provide the full inverse of the sparse matrix – you only get a very specific subset!

5.

Option --hchol carries out the Cholesky factorisation of a sparse matrix only. Calculations are similar to the first step invoked by option --invspa but the matrix is processed in the original order, i.e. no reordering to minimize fill-in is carried out.
Output: The Cholesky factor is written out to filename.chlsky with one row per non-zero element, containing row number, column number and the element of the matrix (space separated).

6.

Option --hdiag provides another variant of --invspa to read a *.gin file (allowing for the header line) and write out the diagonal elements of the inverse only. This is provided with the calculation of the diagonal elements of H from H⁻¹ for single step analyses in mind. Note that it treats the whole of H⁻¹ as sparse – if a high proportion of animals are genotyped, inversion treating the matrix as dense may be more appropriate.
Output: Diagonal elements are written out to filename.hdiags with each row containing row number, original ID (transferred from corresponding .codes file, if available) and the diagonal element of the matrix (space separated).

7.

Option --invspappxn₁,n₂,n₃ provides an approximation of the diagonal elements only of the inverse of a large, sparse matrix. It is intended for matrices too large to invert directly. This uses a Gibbs sampler [20] with n₁ samples in total, a burn-in phase comprising n₂ samples and a thinning rate of n₃. Default values are 55,000, 5,000 and 50. Omitting either n₁,n₂,n₃, n₂,n₃ or n₃ is permissible and will accept the respective defaults.
Output: Diagonal elements are written out to filename.diainv with each row containing the row number and diagonal element of the matrix (space separated).

For Linux versions of WOMBAT (compiled using ifort), the LAPACK routines used are loaded from the Intel MKL library which is parallelised and highly optimised.

Output of inverse: The inverse is written out to filename.inv with one row per non-zero element, containing row number, column number and the element of the matrix (space separated).

EXAMPLE:

wombat -v ---amd ---invspa matrix.dat

specifies sparse matrix inversion of the matrix stored in matrix.dat, using approximate minimum degree ordering and requesting ‘verbose’ output. The output file containing the inverse is matrix.dat.inv.

5.2.11 Calculation of genomic relationships or H⁻¹

Run option --hinv invokes a pre-analysis step to calculate the (inverse of) the joint relationship matrix between genotyped and  non-genotyped animals in a single-step analysis. There are numerous options to select output, alignment and subsets of calculations.

Detailed documentation is not included in this manual, but a pdf file describing the various options is available in Example 20 (chapter 9).

5.2.12 Quadratic approximation

For some types of analyses a parameter ( e.g. the dilution factor for “social” genetic effects) is fixed at a given value and variance components are estimated for this value. To estimate the parameter then requires multiple runs for different values, each generating a point of the profile likelihood for this parameter. Provided an initial triplet of points can be established, comprised of a ‘middle’ parameter value with the highest likelihood bracketed by points with lower likelihoods, a quadratic approximation of the curve can be used to locate its maximum. WOMBAT collects the pairs of parameter values and corresponding likelihoods in a file with the standard name LogL4Quapprox.dat. Specifying the run option --quapp invokes an attempt at a quadratic approximation, utilizing the information from this file. If successful (i.e. if a suitable triplet of points can be found), the estimate of the parameter value at the maximum of the parabola is given, together with its approximate standard error.

5.2.13 Analyses of subsets of traits

For multivariate problems involving more than a few traits, it is often desirable to carry out analyses considering subsets of traits, for example, to obtain good starting values or to trace problems in higher-variate analyses back to particular constellations of traits.

WOMBAT provides the run time option --subsetn to make such preliminary analyses less tedious. It will cause WOMBAT to read the parameter file for a multivariate analysis involving q traits and write out the parameter files for all q uni- (n = 1) or q(q − 1)∕2 bi-variate (n = 2) analyses possible. For n > 2, the program will prompt for the running numbers of the n traits to be considered and write out a single parameter file. These subset analyses are based on the data (and pedigree) file for the ‘full’ multivariate model, using the facility for automatic renumbering of traits and subset selection, i.e. no edits of data files are necessary.

N.B.: This option does not carry through any information from a SPECIAL block, and is not implemented for random regression models or analyses involving correlated random effects or permanent environmental effects fitted as part of the residuals.

For n = 2 (bivariate), WOMBAT also writes out a file containing a bash shell script, run_bivar.sh, to allow execution of all analyses with a single command.

On analysis, encountering the syntax for trait renumbering (see subsection 4.10.7) causes WOMBAT to write out an additional file with the estimates for the partial analysis (Standard name EstimSubsetn + … + m.dat with n to m the trait numbers in the partial analysis, e.g. EstimSubset2+7.dat; see subsection 7.2.6). In addition, the name of this output file is added to a file called SubSetsList.

HINT: WOMBAT will add a line to SubSetsList on each run – this may cause redundant entries. Inspect & edit if necessary before proceeding to combining estimates !

5.2.14 Pooling estimates of covariance components from part analyses

Combining estimates of covariances from analyses involving different subsets of traits is a regular task. This may be a preliminary step to a higher-dimensional multivariate analysis to obtain ‘good’ starting values. Alternatively, analyses for different subsets of traits may involve different data sets – selected to maximise the amount of information available to estimate specific covariances – and we may simply want to obtain pooled, possibly weighted estimates of the covariance matrices for all traits which utilise results from all partial analyses and are within the parameter space.

WOMBAT provides two procedures to combine estimates of covariance components from different analyses or modify existing matrices. These options are not available for analyses involving random regression models, correlated random effects or permanent environmental effects fitted as part of the residuals.

5.2.14.1 Iterative summing of expanded part matrices

Option --itsum selects a run to combine estimates from partial analyses, using the ’iterative summing of expanded part matrices’ approach of Mäntysaari [27] (see also Koivula et al. [23]), modified to allow for differential weighing of individual analyses. For this run, the parameter file required is equal to that for a full multivariate analysis for all traits. Further, a file SubSetsList is assumed to exist and list the names of all the files containing the results from analyses of subsets, and, optionally, the weightings to be applied (see subsection 7.3.9). Pooled covariance matrices are written to a file name PDMatrix.dat as well as a file named PDBestPoint (see subsection 7.2.7).

HINT: Use of --itsum is not limited to combining bi-variate analyses or the use of files with standard names (EstimSubsetn+…+m.dat), but all input files must have the form as generated by WOMBAT.

To use --itsum to combine estimates from analyses involving different data sets, be sure to a) number the traits in individual analyses appropriately (i.e. 1,…,q with q the total number of traits, not the number of traits in a partial analysis), and b) to use the syntax described in subsection 4.10.7 to renumber traits – this will ‘switch on’ the output of subset results files.

Copy PDBestPoint to BestPoint and run WOMBAT with option --best to obtain a listing with values of correlations and variance ratios for the pooled results.

5.2.14.2 Pooling using a (penalized) maximum likelihood approach

Option --pool is similar to --itsum but a) employs a maximum likelihood approach, as described by Meyer [34] , b) facilitates pooling of covariance matrices for all sources of variation simultaneously, and c) allows for penalties to be imposed aimed at ‘improving’ estimates by reducing sampling variation.

Input is as for --itsum and pooled estimates are written to files named PoolEstimates.out (summary file) and PoolBestPoint.

HINT: Details and a file with further documentation are available with Example 15.

5.2.15 Miscellaneous

Option --expiry will print the expiry date for your copy of WOMBAT to the screen.

Option --limits can be used to  find at the upper limits imposed on analyses feasible in WOMBAT, as ‘hard-coded’ in the program. N.B. Sometimes these are larger than your computing environment (memory available) allows.

Option --times causes WOMBAT to  print out values for the CPU time used in intermediate steps.

Option --wide will generate formatted output  files which are wider than 80 columns.

Run option --help causes a list of available run options to be  printed to screen.

5.3 Advanced run options

There are default values for most of these codes, which are adequate for most simpler analyses, and analyses of small to moderately sized data sets. Hence these options are predominantly of interest to users which want to fit complicated models or analyse large data sets, and thus need to tweak some of the options for ordering of equations, parameterisations, maximisation strategies and convergence criteria to get the best possible performance from WOMBAT.

In contrast to the basic options above, several of these options can be followed by one or more, optional numerical arguments. If such arguments are given, they must follow immediately (no spaces), and multiple options need to be separated by comma(s). Arguments must be given sequentially, i.e. if we want to set argument k, all preceding arguments (1,…,k − 1) are required as well. If no arguments are given, the corresponding variables are set to their default values (see Table A.1).

5.3.1 Ordering strategies

The order in which the equations in the mixed model are processed can have a dramatic impact on the computational requirements of REML analyses. Hence, WOMBAT attempts to reorder the equations, selecting a default ordering strategy based on the number of equations; see section A.1 for details. This can often be improved upon by selecting the strategy used explicitly.

Option --mmd specifies ordering using the multiple minimum degree procedure.

Option --amd selects an approximate  minimum degree ordering [2].

Option --metis selects an ordering using a multilevel nested dissection procedure. Up to four optional arguments modifying the behaviour of this subroutine can be specified.

1.

The number of graph separators to be considered, which can be between 0 and 20. As a rule, the higher this number, the better the quality of ordering tends to be. However, the time required for ordering increases substantially with this number, and in some cases intermediate values (between 3 and 9) have been found to work best. The manual for MeTis recommends values up to 5. However, Meyer [31] found values as high as 12 or 14 to be advantageous for large analyses. By default, WOMBAT uses a value of 5 for analyses with more than 50000 and up to 200000 equations, and a value of 10 for larger models.

2.

A factor which tells MeTis which rows in the matrix are to be treated as dense. For instance, a value of 150 means all rows with more than 15% (factor divided by 10) elements than average. The default used by WOMBAT is 0.

3.

An option to select the edge matching strategy employed by MeTis. Valid values are 1 for random, 2 for heavy and 3 for standard edge matching. WOMBAT uses a default of 1.

4.

An option to select whether MeTis uses one- or two-sided node refinement. WOMBAT uses a default of 1 which, paradoxically (to agree with the MeTis convention), invokes two-sided refinement. This is deemed to produce somewhat better orderings. An option of 2 here selects one-sided refinement, which can be faster.

In addition, the option --metisall is available. This causes WOMBAT to cycle through the number of graph separators from 5 to 16, considering density factors of 0 and 200, as well as random and standard edge matching (48 combinations). Warning : This can be quite time-consuming !

To obtain the ordering, WOMBAT assigns a large matrix to store information on the non-zero elements  in the mixed model matrix. The default size chosen is based on the number of equations and traits analysed, and is meant to be generous enough to be sufficient for a wide range of cases. In some instances, however, this can result in WOMBAT trying to allocate arrays which exceed the amount of RAM available and thus failing to run while, in reality, much less space is required for the analysis concerned. In other cases, the default guess is simply not sufficient. To solve these problems, the run time option --choozhz is supplied, which allows the user to set this number. If used as shown, it causes WOMBAT to pause, write out the default value, and read in (input from the terminal!) the new value. Such interactive input may be inconvenient in some scenarios. Hence, alternatively, the new value may be specified as part of the run time option, immediately after (no spaces!) --choozhz. For example, --choozhz11000000 sets the value no 11 million. If --choozhz0 is given, WOMBAT will replace the 0 with the current, hard-coded maximum value for the array size (use run option --limit to find out what it is); this is convenient if something larger than the default is required and plenty of RAM is available.

5.3.2 REML algorithms

WOMBAT can be told to use a particular algorithm to locate the maximum of the likelihood function. In the following, let n (must be an Integer number) denote the maximum number of iterates to be carried out by an algorithm, and let C denote the convergence criterion to be used. Unless stated otherwise, C represents the threshold for changes in log likelihood between subsequent iterates, i.e. convergence is assumed to be reached if this is less than C.

Option --aireml specifies a ‘straight’ AI algorithm. It can be followed by values for n and C. Valid forms are --aireml, --airemln and --airemln,C but not --airemlC.

EXAMPLE:

--aireml30,0.001

limits the number of iterates carried out to 30, and stops estimation when the change in log likelihood is less than 10⁻³.

The computationally most expensive part of an AI iterate is the calculation of first derivatives of the log likelihood. Two methods are available: --derinv selects calculation based on the sparse inverse of the coefficient matrix (default) – this tends to be advantageous when multiple threads are available for parallel execution. --deraut specifies evaluation using automatic differentiation of the Cholesky factor of the coefficient matrix; this was the previous default and is maintained for compatibility.

By default, the AI algorithms enforces an increase in log likelihood in each iterate. This can be switched off using the option --nostrict. This can improve convergence in some cases.

In this case, a modification of the AI matrix can result in better performance of the AI algorithm. Four procedures have been implemented in WOMBAT. These are selected by specifying --modaimn, with n = 1,2,3,4 selecting modification by setting all eigenvalues to a minimum value (n = 1), by adding a diagonal matrix (n = 2, the default), or through a modified (n = 3) [45, 46] or partial (n = 4) [15] Cholesky decomposition of the AI matrix; see section A.5 for details.

Option --emalg selects estimation using a standard EM algorithm. As for --aireml, this option can be followed by n or n,C to specify the maximum number of iterates allowed and the convergence criterion.

Option --pxem then selects estimation using the PX-EM algorithm. As above, this can be followed by n or n,C.

Option --pxai specifies a hybrid algorithm consisting of a few initial rounds of PX-EM, followed by AI REML. This is the default for full rank estimation (unless -c is specified without any explicit definition of the algorithm to be used). This option can have up to three sequential arguments, m denoting the number of PX-EM iterates, and n and C as above.

EXAMPLE:

---pxai6,50,0.001

defines 6 PX-EM iterates, followed by up to 50 AI iterates and a convergence criterion of 10⁻³ for the log likelihood.

Option --emai is like --pxai except that the initial iterates are carried out using the standard EM algorithm. This is the default for reduced rank estimation. Optional arguments are as for --pxai.

Option --cycle, followed by an optional argument n, is used in conjunction with --pxai or --emai. If given, it will repeat the number of (PX-)EM and AI iterates specified by these options n times. If n is omitted, the number of cycles is set to 100. This option is useful for reduced rank analyses, where cycling algorithms appears to be beneficial.

Finally, WOMBAT incorporates two choices for a derivative-free algorithm. While little used, these can be usefully to check for convergence in ‘tough’ problems. Option --simplex selects the simplex or polytope algorithm due to Nelder and Mead [41], as previously implemented in DfReml. This option can have three arguments, n and C as above (but with the convergence criterion C describing the maximum variance among the log likelihood values in the polytope allowed at convergence), and S the initial step size. Option --powell invokes maximisation using Powell [43]’s method of conjugate directions, again ported from DfReml and with optional parameters n, C (as for --aireml) and S.

Option --force can be used to enforce a check for convergence after AI iterates have terminated. It invokes what is essentially a continuation run with --powell followed potentially by more AI iterates. This is done automatically for a small number of selected analyses, but not part of the normal procedure as it can be computationally demanding.

Not really a REML algorithm: Option --like1 selects a single likelihood evaluation for the current starting values, or, if combined with -c, the currently best point as contained in BestPoint. When combined with -v, the individual components of log ℒ are printed to the standard output (screen).

Similarly, option --valid does nothing more than calculate likelihood values. It is provided to aid in the use of cross-validation to estimate the tuning factor for penalised estimation. This involves obtaining estimates for a range of tuning factors for a ‘training’ data set. In doing so, WOMBAT creates a file PenBestPoints.dat (see subsection 7.2.11). The corresponding likelihood values in the ‘validation’ data set are then easily obtained with the --valid run time option: for each set of estimates in PenBestPoints.dat, WOMBAT calculates the likelihood and writes it together with the tuning factor to a file ValidateLogLike.dat.

5.3.3 Parameterisation

n

By default, WOMBAT reparameterises to the elements of the Cholesky factor of the covariance matrices to be estimated (except for full rank analyses using the PX-EM or EM algorithm, which estimate the covariance components directly).

As a rule, the Cholesky factorisation is carried out pivoting on the largest diagonal element. This can be switched off with the --noreord option.

Reparameterisation to the elements of the Cholesky factors removes constraints on the parameter space. Strictly speaking, however, it leaves the constraint that the diagonal elements should be non-negative. This can be removed by transforming the diagonal elements to logarithmic scale. This is selected using the option --logdia, which applies the transformation to all covariance matrices. Conversely, the option --nologd prevents WOMBAT from carrying out this transformation. If neither option is set (the default) and WOMBAT encounters small diagonal elements in any covariance matrices to be estimated during iterations, it will switch to taking logarithmic values of the diagonal elements for the respective matrices only.

5.3.4 Matrix storage mode

By default, WOMBAT assumes that the mixed model equations are sparse and carries out most computations using sparse matrix storage. In certain cases, this may be inappropriate and lead to substantial overheads – a typical example is the case where a random effect has a user-defined covariance matrix which is dense (e.g. a genomic relationship matrix). For this case, the option --dense is provided – specifying this option causes the one triangle of coefficient matrix to be stored assuming all elements are non-zero. Similarly, all calculations are carried out without checking for zero elements - this is done loading subroutines from the BLAS [9] and LAPACK [3] libraries, using packed storage of the coefficient matrix. For n equations, the latter requires n(n + 1)∕2 elements to be stored. With integer*4 addressing,  this sets an upper limit of n = 65,535. For larger problems, the option --dense2 is available which uses storage in two-dimensional arrays instead. Note that use of --dense2 may require substantial amounts of RAM and that it is not yet implemented for all algorithms (e.g. not for PX(EM)).

5.3.5 Sparse matrix factorisation, auto-differentiation and inversion

By default, WOMBAT now employs a ‘super-nodal’ approach in the sparse matrix manipulations. This involves gathering dense sub-blocks and processing these using BLAS [9] and LAPACK [3] library routines. While this readily allows processing  utilising multiple cores in these steps, the library routines require matrices to be strictly positive definite. Occasionally, they may fail where the previous procedure used may have circumvented the problem. The run time option --old is available to switch back to the latter.

5.3.6 Other

5.3.6.1 Pedigree pruning

By default, WOMBAT ‘prunes’ the pedigree information supplied, i.e. eliminates any uninformative individuals (without records and links to only one other individual), before calculating the inverse of the  numerator relationship matrix. Exceptions are prediction runs (option --blup) and models which fit correlated additive genetic effects. In some instances, however, it may be desirable not to ‘prune’ pedigrees. Pruning can be switched off through the run time option --noprune.

5.3.6.2 Overriding defaults for pedigree reduction

For sire model analyses, the default is not to carry out any pruning or pedigree reduction. The option --redped is provided to switch on pedigree reduction for sire model analyses, i.e. elimination of all individuals without a link to an individual in the data from  the pedigree file.

Occasionally – for instance to compute breeding values for progeny of animals in the data – it is desirable not to reduce the pedigree prior to analysis. Providing the  option --norped eliminates this step and switches off ‘pruning’ at the same time. Option --redped is invoked automatically for runs with options --solvit or --s1step.

5.3.6.3 No programmed pause please

Generally, WOMBAT does not require any interactive input. An exception is a ‘pause’ after a warning message on a constellation of options which is new or can be problematic. These programmed pauses can be switched off using the run time option --batch.

5.3.6.4 Centering data

By default, WOMBAT ‘centers’ data, i.e. subtracts the trait-specific mean from each observation. For some (rare) tasks this is undesirable. Hence, the option --nocenter is provided which eliminates this step.

5.3.6.5 Choosing the algorithm to calculate inbreeding

By default, WOMBAT uses the algorithm of Tier [47] to compute inbreeding coefficients before setting up the inverse of the numerator relationship matrix between individuals. This algorithm is implemented with a limit of 23 generations. Occasionally, this results in insufficient space, both for pedigrees with mode generations and very large pedigrees. The option --maxgen, followed immediately by two digits giving the new value (e.g. maxgen25) allows this limit to be increased, up to a maximum value of 30. Note that this may result in WOMBAT trying to assign some very large arrays, and may slow calculations.

Alternative algorithms, due to Quaas [44] and Meuwissen and Luo [28], which do not have this limitation but are slower and use very little memory, can be selected using the run time option --quaas and --meuw, respectively.

5.3.6.6 Using only the first n observations

Sometimes, it is useful to carry out a preliminary analysis considering a small(ish) subset of the data only. For instance, we may want to test that we have fitted the model correctly or we may want to estimate suitable starting values. A quick way to do this is available through the option --nonlyN, where N is an integer number specifying the number of records to be used (no space between nonly and N). For example, --nonly1000 selects an analysis utilising the first 1000 records only.

5.3.6.7 Omitting the preliminary LSQ analysis

One of the preliminary steps carried out by WOMBAT is a simple least-squares analysis of the fixed part of the model – this aims at identifying any dependencies amongst fixed effects levels which might not have been correctly specified. However, if there are many cross-classified effects or covariables, this can be very slow, especially for large data sets (code for this step is very old and probably the weakest part of WOMBAT and not as highly  optimised as the main part of the calculations). Option --nolsq is available to skip this step and is recommended for large analyses with many cross-classified effects – but it does rely on all rank deficiencies in the part of the coefficient matrix for fixed effects being identified: the run is likely to fail with an appropriate error message otherwise.

A scenario where this option is useful is when we fit different random effects for a given data set where the fixed part of the model remains constant – typically, it is then sufficient to carry out the least-squares step only for the first run. Alternatively, it is beneficial for simulations studies where the data is re-sampled, so that data and pedigree structure are the same for all replicates.

5.3.6.8 Omit writing out solutions

Run option --nosolut will skip writing out  files with fixed and random effects solutions at the end of an estimation run. This is useful for simulation studies where these are not of interest.

5.3.6.9 Headers for solution files

Run option --nohead will omit writing out  the header line for files with solutions for random effects. This is useful when these files are to be processed further using other software.

5.3.6.10 Raw mean

Run option --noraw will omit writing out  the raw means for solutions of fixed or random effects. This is useful for large problems and when these are not of interest.

5.3.6.11 Selecting the number of threads for execution

The Linux executables compiled using ifort will carry out calculations using multiple threads where available. The default is to use half of the maximum number available (as each core typically is counted as two threads, this is the number of cores). A different number can be selected via the run option --threadsn where n is an integer value giving the number of threads to be utilised.

NB: An alternative is to set the maximum number using the environment variable OMP_NUM_THREADS – WOMBAT will use half of that number (default) or up to that number if --threads is specified.

5.4 Parameter file name

If a parameter file other than wombat.par is to be used, this can be given as the last entry in the command line. The extension .par can be omitted.

EXAMPLE:

wombat  weights

specifies a run where WOMBAT expects the parameter file weights.par.