WOMBAT offers a substantial number of so-called advanced run time options. While simple analyses tend to work well without them, using these options can greatly enhance reliability and efficiency for analyses involving more than a few traits (or RR coefficents) or multiple random effects. This is especially important for large data sets where each additional iterate or likelihood evaluation can require substantial time !
Choice of parameterisation can have a big impact on the convergence rates of REML analyses.
By default, using the average information (AI) or reduced rank Expectation Maximisation (EM or PX-EM) algorithm, WOMBAT estimates the elements of the Cholesky factors of the covariance matrices specified rather than the covariance components themselves. This is sufficient to avoid problems with estimates `out of the parameter space' for a substantial number of cases. However, this does not provide a completely unconstrained parameterisation: For the corresponding covariance matrices to be positive definite, the diagonal elements of the Cholesky decomposition must be greater than zero.
Hence, for analyses where some of the diagonal elements are expected to become small, it is advisable to estimate their log values, thus removing any restriction. This parameterisation is enforced with the run option −−logdia.
A typical example where this is appropriate is the analysis of highly correlated traits, such as weights at different ages, in particular if maternal effects are fitted.
WOMBAT does not invoke this parameterisation to begin with, as convergence tends to be better without it if small diagonals are not an issue. While WOMBAT does some damage control – automatically switching to log diagonals if very small diagonals are encountered – performance is generally better if this is specified explicitly to begin with for analyses involving highly correlated random effects. This holds especially for full rank analyses. For reduced rank analyses, omitting all principal components with small eigenvalues is generally sufficient to eliminate problems with diagonal elements in the Cholesky factor close to zero.
Option −−nologd prevents automatic switching to log diagonals during the course of AI iterates.
When determining the Cholesky decomposition of the covariance matrices to be estimated, WOMBAT pivots on the largest diagonal, i.e. considers traits in descending order of variance, given all traits processed so far. The mixed model equations are set up and the re-ordering of equations is determined for this sequence, derived from the matrix of starting values.
It would be easiest to maintain this order for a particular analysis. If not, we may need to repeat part of the set-up steps, i.e. finding a permutation of the original equation numbers to minimise `fill-in' and operations per iterate and determining the sparse matrix structure of the mixed model matrix (symbolic factorisation).
However, in some instances it is advantageous – for speed and reliability of convergence – to allow the order to be changed. Most often, this applies to multivariate analyses with `bad' starting values or reduced rank analyses. Bear in mind though that here is a clear trade-off: for large analyses, the additional time required to repeat the set-up and reordering steps can be substantial.
−−norecon (for `no reconstruction' of the mixed model matrix) specifies a run where no re-ordering is allowed. This is the default for full rank, multivariate analyses.−−reconstr does the opposite, allowing for potential reordering. This is the default for all reduced rank analyses. REML analyses in WOMBAT are carried out using an iterative solution scheme. Iterations continue until a predefined maximum number of iterates is reached or when changes in estimates or likelihood between subsequent iterates become sufficiently small.
The default maximum number of iterates allowed is fairly generous and default convergence criteria used by WOMBAT are quite stringent.
Hence, without additional control often more iterates are performed than are strictly necessary or desired. This holds especially for `prelimary' analyses (where we don't really need results to be all that accurate), or large problems (where additional iterates represent a large computational burden).
Control of the number of iterates performed is linked to the run time options to select the algorithm used to find the maximum of the likelihood function.
The default number of iterates used depends on the number of parameters to be estimated. Additional iterates are allowed for reduced rank analyses. To override the default, the maximum number of iterates needs to be specified as the first optional part of a run time option selecting a maximisation algorithm. For instance:
−−aireml44 selects the AI-REML algorithm and allows a maximum of 44 iterates−−emalg1000 specifies maximisation via the EM algorithm with a maximum of 1000 iterates−−pxai12,20 specifies maximisation using 12 iterates of the PX-EM algorithm, followed by up to 20 iterates using the AI step.By default, WOMBAT continues iterations until the chance in log likelihood between subsequent iterates drops below 0.0005 for the AI algorithm and below 0.00001 for the (PX-)EM algorithm. If it is desired to use a different limit, this can be specified together with the option for the maximisation algorithm and the maximum number of iterates. It needs to be given after the number of iterates, separated from it by a comma (no spaces). For instance
−−aireml18,0.001 specifies a run using the AI REML algorithm with a maximum of 18 iterates, stopping when a change in log likelihood between iterates of less than 0.001 has been encountered.−−pxai3,11,0.005 chooses 3 iterates of the PX-EM algorithm, followed by up to 11 AI iterates, stopping the change in log L becomes less than a value of 0.005.WOMBAT offers a choice of maximisation algorithms. Again, the defaults provided are suitable for a wide rannge of problems, but for specific cases it is desirable to control the algorithm(s) used.
The main work-horse is the `average information' (AI) algorithm. This is the default for analyses involving up to 18 parameters to be estimated (i.e. most uni- and bivariate analyses), and any continuation runs (run option −c).
For other analyses, use of the AI algorithm can be enforced with the command line option −−aireml.
The AI algorithm represents a modified Newton-Raphson type algorithm to locate the maximum of the likelihood function. It is generally fast to converge. However, problems have been encountered, in particular for
While WOMBAT attempt to exert strict control over the AI steps, it may be advantageous to carry out a few iterates using an expectation maximisation algorithm for such cases.
Expectation Maximisation (EM) type algorithms for mixed model analyses are renown for their reliability but also their slow convergence. A particular variant is the 'parameter expanded' EM (PX-EM) algorithm. This has been shown to converge considerably faster than the standard EM algorithm.
For the first few iterates in a particular analysis, the PX-EM algorithm often yields larger steps towards the maximum of the log likelihood than the AI algorithm. Hence it appears sensible to carry out the initial iterates using the PX-EM algorithm and to switch to the AI algorithms subsequently. However:
Use of the PX-EM algorithm is enforced by the command line option −−pxem.
A combination of a PX-EM steps followed by b AI iterates is selected by the option −−pxaia,b, e.g. −−pxai10,50.
For analyses with more than 18 parameters to be estimated, the default implemented in WOMBAT is that the first 3 iterates are PX-EM iterates.
This can be modified by option −−good (referring to 'good' starting values) which reduces he number of PX-EM steps to 1. Conversely, option −−bad increases the PX-EM steps to up to 8, and doubles the default number of subsequent AI steps (the switch to AI may occur after less than 8 PX-EM iterates if the change in likelihood between iterates becomes less than 2).
WOMBAT offers two derivative-free maximisation procedures:
Both are a legacy from DFREML and of relatively little practical use. Ocacsionally, they come in handy to check for convergence or to help when all the other alternatives fail.