5.2 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 A.1).


Table 5.2: ‘Advanced’ run options for WOMBAT
Option

Purpose

Ordering strategies for mixed model matrix (MMM)
--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

Specific REML algorithms
--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 definite

--emalg

Use the standard EM-algorithm

--pxem

Use the PX-EM algorithm

--emai

Use a few rounds of EM followed by AI REML (default)

--pxai

Use a few rounds of PX followed by AI REML

--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)

--like1

Carry out a single likelihood evaluation only

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

Miscellaneous
--noprune

Do not prune pedigrees

--batch

Do not pause after warning messages

--quaas

Use alternative algorithm [24] to set up NRM inverse

--redped

Switch on reduction of pedigree file for a sire model.

5.2.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 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 [1].

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 [17] 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 or accessible (some of the 32-bit versions are restricted to 2GB) and thus failing to run while, in reality, much less space is required for the analysis concerned. Hence, the run time option --choozhz is supplied: This cases WOMBAT to pause, write out the default settings, and read in (from the terminal) an alternative value specified by the user.

5.2.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−3  .

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  ) [2526] or partial (n = 4  ) [8] Cholesky decomposition of the AI matrix; see 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: --pxem6,50,0.001
defines 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 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 three 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 [20], 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 [22]’s method of conjugate directions, again ported from DFREML and with optional parameters n  , C  (as for --aireml) and S  .

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

5.2.3 Parameterisation

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.2.4 Other

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.

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.

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.

By default, WOMBAT uses the algorithm of Tier [27] to compute inbreeding coefficients before setting up the inverse of the numerator relationship matrix between individuals. This algorithm is implemented with a hard-coded limit of 23 generations. An alternative algorithm, due to Quaas [24], which does not have this limitation but is a lot slower, can be selected using the run time option --quaas.