Run options

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 eﬀects of options given ﬁrst. 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 ﬁle with the default name RunOptions in the current working
directory.

This ﬁle must contain a single option per line. If such ﬁle exists, it is read before any command line options are parsed. There are no checks for options speciﬁed 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 ﬁrst line of the parameter ﬁle, after a code of RUNOP (see section 4.3). Any such options are read before options speciﬁed 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.

Option | Purpose [Section] |

-c | Specify a continuation run [5.2.1 ] |

-v | Specify verbose screen output [5.2.2 ] |

-d | Specify very detailed screen output [5.2.2 ] |

-t | Specify terse screen output [5.2.2 ] |

Default REML algorithms
| |

--good | Tell that you have good starting values [5.2.4 ] |

--bad | Tell that you have bad starting values [5.2.4 ] |

Non-estimation runs
| |

--setup | Select a run performing the set-up steps only [5.2.3 ] |

--best | Select a run printing out estimates for the currently best point only [5.2.6 ] |

--blup | Select a prediction (BLUP) run; direct solution [5.2.7 ] |

--solvit | Select a prediction (BLUP) run; solve iteratively [5.2.7 ] |

--mmeout | Like --solvit, but write MME out to ﬁle. [5.2.7 ] |

--s1step | Like --solvit, but for single step model. [5.2.7 ] |

--s2step | Like --s1step, but using iteration on data. [5.2.7 ] |

--hinv | Select calculations of related to H |

--snap | Carry out a GWAS type analysis for multiple SNPs [5.2.7 ] |

--simul | Select a run simulating data only [5.2.8 ] |

--sample | Select a run sampling estimates of covariance matrices [4.17.8 ] |

--subset | Write out parameter ﬁles for analyses of subsets of traits [5.2.13 ] |

--itsum | Iterative summation of partial covariance matrices [5.2.14.1 ] |

--pool | Pooling covariance components by penalized maximum likelihood [5.2.14.2 ] |

--invert | Invert a dense symmetric p.d. matrix; use LAPACK routines [5.2.10 ] |

--inveig | Invert a dense symmetric matrix; use eigen decomposition [5.2.10 ] |

--invlap | Invert a dense symmetric p.d. matrix; use LAPACK routines [5.2.10 ] |

--invspa | Invert a sparse symmetric matrix[5.2.10 ] |

--hdiag | Invert a sparse symmetric matrix; write out diagonals only[5.2.10 ] |

--hchol | Carry out Cholesky factorisation of a sparse symmetric matrix [5.2.10 ] |

--quapp | Perform quadratic approximation of likelihood [5.2.12 ] |

Matrix storage
| |

--dense | Specify dense mixed model equations [5.3.4 ] |

Miscellaneous
| |

--expiry | Print out the expiry date for the program (to screen) [5.2.15 ] |

--limits | Print out the hard-coded program limits (to screen) [5.2.15 ] |

--times | Print out various intermediate times for a run [5.2.15 ] |

--wide | Write ‘wide’ output ﬁles [5.2.15 ] |

--nohead | Omit header line for RnSoln_*.dat ﬁles [5.2.15 ] |

Numerical settings
| |

--zero | Change value of operational zero [5.2.5.1 ] |

--pivot | Change value for smallest pivot in Cholesky factor (covariance matrix) [5.2.5.2 ] |

--fudge | Add small constant to diagonals of coeﬃcient matrix [5.2.5.3 ] |

In some instances, REML estimation continuing from the ‘best’ estimates so far is
necessary. This is invoked with the -c option.

If speciﬁed, WOMBAT will attempt to read its ‘starting’ values from the ﬁle
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 ﬁle instead. The -c option selects estimation using the AI algorithm,
unless another procedure is selected explicitly.

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 ﬁle
- -d
- selects detailed output, useful for debugging

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 ﬁle SumModel.out (see subsection 7.1.2) is as expected, i.e. that WOMBAT is ﬁtting the correct model and has read the data ﬁle correctly. This is speciﬁed using --setup. This option also invokes verbose screen output.

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 conﬁdent, 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.

WOMBAT uses a small positive value as operational zero to reduce numerical
errors. The default value is 10^{−8}. This can be changed specifying the run option
--zeron where n is a single-digit integer. This option redeﬁnes the operational zero
to be 10^{−n}.

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.

In some instances a run can fail because of numerical issues in the Cholesky
factorisation of the coeﬃcient matrix of the MME. Most commonly this involves
LAPACK routine DPOTRF. Reasons may include a build up of rounding errors
or a coeﬃcient matrix which is not safely positive deﬁnite to begin with.
A possible work around then is to add a small, positive constant to the
diagonal elements of the coeﬃcient matrix. This can be requested by adding

--fudge

to the command line (or ﬁrst line in the parameter ﬁle). The default is to add 10^{−8}.
This can be changed to 10^{−n} by specifying --fudgen instead with n an appropriate
integer value.

WOMBAT writes out the currently best values of estimates of covariance components to the ﬁle BestPoint whenever the likelihood is increased. The option --best causes WOMBAT to read this ﬁle and write out the matrices of covariances and corresponding correlations in more readily legible format to the ﬁle 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 ﬁles used.

WOMBAT can be used for a simple BLUP run, using the ‘starting’ values given in the parameter ﬁles as assumed values for the true covariances. In this mode, no pruning of pedigrees is carried out. If -c is speciﬁed in addition, WOMBAT will try to read the values from the ﬁle 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 coeﬃcient 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 eﬀects) for random eﬀects to be calculated from the diagonal elements of the direct inverse.
- Run option --solvit speciﬁes 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 eﬀects
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
^{−8}. 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^{−6}.

This option increases the limit on the maximum number of eﬀects (equations) which can be ﬁtted 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 ﬁle containing all non-zero elements in the lower triangle of the coeﬃcient matrix (default name MMECoeffMatrix.dat), and the ﬁle MMEEqNos+Solns.dat which speciﬁes the mapping of equation numbers to eﬀects ﬁtted, 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 diﬀerence 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
^{−1}) 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 ﬁle.Change: The default for the PCG algorithm has been changed to the so-called SSOR preconditioner. There are now ﬁve 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 eﬀects for all traits in a multivariate analysis. However, if genetic groups are ﬁtted as an additional random eﬀect, 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 eﬀects other than genetic groups (equivalent to B for univariate nalyses).

A B C D E ’diagonals’ ’diagblock’ ’genoblock’ ’diaggenos’ ’diagggrps’ genotyped animals

diag ^{a}blox full full diag non-geno. anim.s

diag blox blox diag diag genetic groups

diag full full full full other eﬀects

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 ﬁxed variance components (i.e. a BLUP run)– ﬁtting a mixed model with a single SNP eﬀect (ﬁtted as a linear covariable) but processing multiple SNPs with reference allele counts read in from a separate ﬁle (default name QTLAllels.dat or QTLAllelsR.dat; see subsection 6.6.3). This assumes no missing counts and employs an eﬃcient computing strategy which exploits that only the equation for the SNP changes for diﬀerent SNPs [40]. Estimated SNP eﬀects together with their standard errors are written out to a ﬁle with the default name QTLSolutions.dat (see subsection 7.2.10). This run option must be used in conjunction with a SPECIAL statement which identiﬁes which of the covariables in the model of analysis is to be treated as SNP eﬀects (see subsection 4.17.7).

Specifying --simul causes WOMBAT to sample values for all random eﬀects ﬁtted for the data and pedigree structure given by the respective ﬁles, from a multi-variate normal distribution with a mean of zero and covariance matrix as speciﬁed in the parameter ﬁle. Again -c can be used to acquire population values from the ﬁle BestPoint instead. No ﬁxed eﬀects 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 speciﬁed explicitly in a ﬁle RandomSeeds (see subsubsection 6.6.5.3). Output ﬁle(s) have the standard name(s) SimData001.dat, SimData002.dat, …, SimDatan.dat. These have the same layout as the original data ﬁle, 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).

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 ﬁles 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 eﬀect only can be selected by specifying its name in a SPECIAL statement (see subsection 4.17.8). Samples generated are written to a ﬁle with the standard name CovSamples_ALL.dat or CovSamples_name.dat with name the name of the random eﬀect selected (see subsection 7.2.12). These ﬁles are ‘formatted’ (text) ﬁles that are suitable for further analysis using standard statistical software packages. In addition, a ﬁle 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 eﬀects with diagonal covariance matrices. It tends to work best for estimation runs with the PC option, even when matrices are estimated at full rank.

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 deﬁnite (see below for exceptions) needs to be supplied in a ﬁle, 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 diﬀerent ‘modes’, invoked by diﬀerent run options, mostly beginning with --inv and followed (space separated) by the name of the ﬁle 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 ﬁle 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 diﬀer 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 ﬁrst performs an eigen-value and -vector decomposition
of the matrix, using LAPACK routine DSYEVR. If eigenvalues less
than a speciﬁed value (default: 10
^{−12}) are found, these are set to this minimum, before constructing the inverse from the eigenvectors and -values. Eigenvalues are written out to a separate ﬁle. - 4.
- Option --invspa chooses sparse matrix inversion. This is suitable for
the inversion of large, sparse matrices where selected elements of the
inverse are suﬃcient. 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 speciﬁc subset!

- 5.
- Option --hchol carries out the Cholesky factorisation of a sparse matrix only.
Calculations are similar to the ﬁrst step invoked by option --invspa but the
matrix is processed in the original order, i.e. no reordering to minimize ﬁll-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 ﬁle
(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
^{−1}for single step analyses in mind. Note that it treats the whole of H^{−1}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 ﬁle, if available) and the diagonal element of the matrix (space separated). - 7.
- Option --invspappxn
_{1},n_{2},n_{3}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_{1}samples in total, a burn-in phase comprising n_{2}samples and a thinning rate of n_{3}. Default values are 55,000, 5,000 and 50. Omitting either n_{1},n_{2},n_{3}, n_{2},n_{3}or n_{3}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`

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

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 ﬁle describing the various options is available in Example 20 (chapter 9).

For some types of analyses a parameter ( e.g. the dilution factor for “social” genetic eﬀects) is ﬁxed at a given value and variance components are estimated for this value. To estimate the parameter then requires multiple runs for diﬀerent values, each generating a point of the proﬁle 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 ﬁle with the standard name LogL4Quapprox.dat. Specifying the run option --quapp invokes an attempt at a quadratic approximation, utilizing the information from this ﬁle. 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.

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 ﬁle for a multivariate analysis involving q traits and write out the parameter ﬁles 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 ﬁle. These subset analyses are based on the data (and pedigree) ﬁle for the ‘full’ multivariate model, using the facility for automatic renumbering of traits and subset selection, i.e. no edits of data ﬁles 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 eﬀects or permanent environmental eﬀects ﬁtted as part of the residuals.

For n = 2 (bivariate), WOMBAT also writes out a ﬁle 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 ﬁle 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 ﬁle is added to a ﬁle 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 !

Combining estimates of covariances from analyses involving diﬀerent 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 diﬀerent subsets of traits may involve diﬀerent data sets – selected to maximise the amount of information available to estimate speciﬁc 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 diﬀerent analyses or modify existing matrices. These options are not available for analyses involving random regression models, correlated random eﬀects or permanent environmental eﬀects ﬁtted as part of the residuals.

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]), modiﬁed to allow for diﬀerential weighing of individual analyses. For this run, the parameter ﬁle required is equal to that for a full multivariate analysis for all traits. Further, a ﬁle SubSetsList is assumed to exist and list the names of all the ﬁles 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 ﬁle name PDMatrix.dat as well as a ﬁle named PDBestPoint (see subsection 7.2.7).

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

To use --itsum to combine estimates from analyses involving diﬀerent 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 ﬁles.

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.

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 ﬁles named PoolEstimates.out (summary ﬁle) and PoolBestPoint.

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

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

Option --limits can be used to ﬁnd 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 ﬁles which are wider than 80 columns.

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

Option | Purpose [Section] |

--metis | Use multi-level nested dissection procedure (MeTis) to re-order MMM |

--mmd | Use multiple minimum degree method to re-order MMM |

--amd | Use approximate minimum degree algorithms for re-ordering |

--choozhz | Assign initial size of mixed model matrix interactively |

--aireml | Use the AI algorithm |

--nostrict | Allow the AI algorithm to take steps decreasing the log likelihood |

--modaim | Choose method to modify the AI matrix to ensure it is positive deﬁnite |

--emalg | Use the standard EM-algorithm |

--pxem | Use the PX-EM algorithm |

--emai | Use a few rounds of EM followed by AI REML |

--pxai | Use a few rounds of PX followed by AI REML (default) |

--cycle | Repeat cycles of (PX-)EM and AI iterates |

--simplex | Use the Simplex procedure for derivative-free (DF) maximisation |

--powell | Use Powell’s method of conjugate directions (DF) |

--force | Force derivative-free search after AI-REML steps |

--like1 | Carry out a single likelihood evaluation only |

--valid | Carry out likelihood evaluations for multiple estimates |

--deraud | Calculate ﬁrst derivatives of log L via automatic diﬀerentiation |

Parameterisations
| |

--logdia | Select log transformation of diagonal elements of Cholesky factor |

--nologd | No log transformation of diagonal elements of Cholesky factor |

--reorder | Pivot on largest diagonal elements of covariance matrices during Cholesky factorisation (default) |

--noreord | Carry out Cholesky decomposition of covariance matrices sequentially |

--reconstr | Allow change in order of pivots of Cholesky factors of covariance matrices and reconstruction of mixed model matrix |

--noreconst | Do not allow changes in order of pivots and reconstruction of mixed model matrix |

Option | Purpose [Section] |

Miscellaneous
| |

--nonly | Consider only the ﬁrst N observations [5.3.6.6 ] |

--noprune | Do not prune pedigrees [5.3.6.1 ] |

--nocenter | Do not subtract mean from observations [5.3.6.4 ] |

--nolsq | Skip least-squares analysis for ﬁxed eﬀects [5.3.6.7 ] |

--batch | Do not pause after warning messages [5.3.6.3 ] |

--maxgen | Increase the space available in Tier”s algorithm to set up NRM inverse |

--quaas | Use algorithm of Quaas [44] to set up NRM inverse [5.3.6.5 ] |

--meuw | Use algorithm of Meuwissen and Luo [28] to set up NRM inverse [5.3.6.5 ] |

--redped | Switch on reduction of pedigree ﬁle for a sire model [5.3.6.2 ] |

--norped | Do not carry out pedigree reduction [5.3.6.2 ] |

--threadsn | Select no. of threads to be used [5.3.6.11 ] |

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 ﬁt 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).

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 speciﬁes 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 speciﬁed.

- 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 reﬁnement. WOMBAT uses a default of 1 which, paradoxically (to agree with the MeTis convention), invokes two-sided reﬁnement. This is deemed to produce somewhat better orderings. An option of 2 here selects one-sided reﬁnement, 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 suﬃcient 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 suﬃcient. 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 speciﬁed 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 ﬁnd out what it is); this is convenient if something larger than the default is required and plenty of RAM is available.

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 speciﬁes 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

^{−3}.

The computationally most expensive part of an AI iterate is the calculation of ﬁrst derivatives of the log likelihood. Two methods are available: --derinv selects calculation based on the sparse inverse of the coeﬃcient matrix (default) – this tends to be advantageous when multiple threads are available for parallel execution. --deraut speciﬁes evaluation using automatic diﬀerentiation of the Cholesky factor of the coeﬃcient 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 oﬀ using the option --nostrict. This can improve convergence in some cases.

In this case, a modiﬁcation 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 modiﬁcation by setting all eigenvalues to a minimum value (n = 1), by adding a diagonal matrix (n = 2, the default), or through a modiﬁed (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 speciﬁes 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 speciﬁed without any explicit deﬁnition 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`

deﬁnes 6 PX-EM iterates, followed by up to 50 AI iterates and a convergence criterion of 10

^{−3}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 speciﬁed 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 beneﬁcial.

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 ﬁle 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 ﬁle ValidateLogLike.dat.

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 oﬀ 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.

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 eﬀect has a user-deﬁned 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 coeﬃcient 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 coeﬃcient 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)).

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 deﬁnite. 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.

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 ﬁt correlated additive genetic eﬀects. In some instances, however, it may be desirable not to ‘prune’ pedigrees. Pruning can be switched oﬀ through the run time option --noprune.

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 ﬁle.

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 oﬀ ‘pruning’ at the same time. Option --redped is invoked automatically for runs with options --solvit or --s1step.

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 oﬀ using the run time option --batch.

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

By default, WOMBAT uses the algorithm of Tier [47] to compute inbreeding coeﬃcients 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 insuﬃcient 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.

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 ﬁtted 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 ﬁrst 1000 records only.

One of the preliminary steps carried out by WOMBAT is a simple least-squares analysis of the ﬁxed part of the model – this aims at identifying any dependencies amongst ﬁxed eﬀects levels which might not have been correctly speciﬁed. However, if there are many cross-classiﬁed eﬀects 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-classiﬁed eﬀects – but it does rely on all rank deﬁciencies in the part of the coeﬃcient matrix for ﬁxed eﬀects being identiﬁed: the run is likely to fail with an appropriate error message otherwise.

A scenario where this option is useful is when we ﬁt diﬀerent random eﬀects for a given data set where the ﬁxed part of the model remains constant – typically, it is then suﬃcient to carry out the least-squares step only for the ﬁrst run. Alternatively, it is beneﬁcial for simulations studies where the data is re-sampled, so that data and pedigree structure are the same for all replicates.

Run option --nosolut will skip writing out ﬁles with ﬁxed and random eﬀects solutions at the end of an estimation run. This is useful for simulation studies where these are not of interest.

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

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

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 diﬀerent 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 speciﬁed.

If a parameter ﬁle 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`

speciﬁes a run where WOMBAT expects the parameter ﬁle weights.par.

Option | Section |

--amd | [5.3.1 ] |

--aireml | [5.3.2 ] |

--bad | [5.2.4 ] |

--batch | [5.3.6.3 ] |

--best | [5.2.6 ] |

--blup | [5.2.7 ] |

-c | [5.2.1 ] |

--choozhz | [5.3.1 ] |

--cycle | [5.3.2 ] |

-d | [5.2.2 ] |

--dense | [5.3.4 ] |

--deraud | [5.3.2 ] |

--derinv | [5.3.2 ] |

--emai | [5.3.2 ] |

--emalg | [5.3.2 ] |

--expiry | [5.2.15 ] |

--force | [5.3.2 ] |

--fudge | [5.2.5.3 ] |

--good | [5.2.4 ] |

--help | [5.2.15 ] |

--hinv | [5.2.11 ] |

--hdiag | [5.2.10 ] |

--hchol | [5.2.10 ] |

--invert | [5.2.10 ] |

--inveig | [5.2.10 ] |

--invlap | [5.2.10 ] |

--invspa | [5.2.10 ] |

--itsum | [5.2.14.1 ] |

--limits | [5.2.15 ] |

--like1 | [5.3.2 ] |

--logdia | [5.3.3 ] |

--maxgen | [5.3.6.5 ] |

--metis | [5.3.1 ] |

--mmd | [5.3.1 ] |

--mmeout | [5.2.7 ] |

--meuw | [5.3.6.5 ] |

--modaim | [5.3.2 ] |

--nohead | [5.3.6.10 ] |

Option | Section |

--noraw | [5.3.6.9 ] |

--nostrict | [5.3.2 ] |

--nologd | [5.3.3 ] |

--noreord | [5.3.3 ] |

--nonly | [5.3.6.6 ] |

--noprune | [5.3.6.1 ] |

--nocenter | [5.3.6.4 ] |

--nolsq | [5.3.6.7 ] |

--norped | [5.3.6.2 ] |

--pivot | [5.2.5.2 ] |

--powell | [5.3.2 ] |

--pxai | [5.3.2 ] |

--pxem | [5.3.2 ] |

--pool | [5.2.14.2 ] |

--quapp | [5.2.12 ] |

--quaas | [5.3.6.5 ] |

--reconstr | |

--redped | [5.3.6.2 ] |

--reorder | |

--s1step | [5.2.7 ] |

--s2step | [5.2.7 ] |

--sample | [4.17.8 ] |

--setup | [5.2.3 ] |

--simplex | [5.3.2 ] |

--simul | [5.2.8 ] |

--snap | [5.2.7 ] |

--subset | [5.2.13 ] |

--solvit | [5.2.7 ] |

-t | [5.2.2 ] |

--threads | [5.3.6.11 ] |

--times | [5.2.15 ] |

-v | [5.2.2 ] |

--valid | [5.3.2 ] |

--wide | [5.2.15 ] |

--zero | [5.2.5.1 ] |