====================================================================================
WOMBAT Example 24: Weighting factors for genetic variances due to SNPs & RPG effects
====================================================================================

This example provides an illustration on how to use WOMBAT to estimate
genetic covariances explained by marker effects and due to residual
polygenic (RPG) effects in a breeding value model.

Provided are simulated data for 3 traits (records, marker counts and
pedigrees), *.par files to build H^{-1} (step 1) and to carry out a
REML or BLUP runs (step 2) as well as bash scripts (run.sh) to run all
steps.

NB: In contrast to the other worked examples, only the input files and
bash scripts "run.sh" are supplied!

A: Univariate REML analysis for trait 1, fitting RPG effects
   explicitly, i.e as an additional random effects. The *.par file
   contains instructions to calculate the total genetic variance, the
   total heritability and the weighting factor lambda as well as their
   approximate standard errors.

B: Analogous to A, but carrying out a multivariate analysis for all
   three traits. Note that a paramterisation to PC scores is specified for
   computational efficiency and that the log determinant of the H matrix
   is calculated to allow model comparisons based on the log likelihood.

C: Multivariate REML analysis for all three traits estimating a
   *single* weighting factor lambda. This involves fitting a single
   genetic effect, representing the sum of marker and RPG effects, and
   accounting for lambda `implicitly' by replacing the GRM in the
   construction of H-inverse by the weighted average of the GRM and the
   NRM for genotyped animals.

   Estimation of lambda is carried out by quadratic approximation of its
   profile log likelihood. Evaluation of each point on this curve
   requires two steps, namely building the H-inverse for the chosen value
   of lambda followed by a REML analysis to determine the pertaining log
   likelihood (logL). To begin with, we evaluate three points, with
   lambdas chosen so that the middle value has the highest logL and is
   bracketed by the other two lambdas. This provides the information
   to fit a quadratic function to the curve -- estimating lambda as that
   corresponding to the maximum of the parabola. After logL for this
   value has been obtained, a new triplet of points is selected and
   another quadratic approximation is carried out to update lambda. This
   need to be repeated until convergence, i.e. until changes in lambda
   and logL become sufficiently small -- typically 3 to 4 iterates
   suffice.

   WOMBAT faciliates the quadratic approximation step by collecting the
   points on the profile likehood in a file named "LogL4Quapprox.dat" and
   providing the run option --quapp to carry out the quadratic
   approximation. In addition, the new value of lambda is written to a
   copy of the *.par file, "wmbnew.par", ready to build Hinverse for the
   next iterate.

D: Multivariate BLUP run (iteration on data) fitting a single genetic
   effect but allowing for trait-specific weighting of genomic and RPG
   effects, exploiting a canonical transformation (illustration only; not
   optimised; undocumented feature in WOMBAT).

D1: Multivariate BLUP run (iteration on data) fitting separate RPG
   effects; showing contrast to D in PCG iterates needed.
