====== Example 1 for WOMBAT ======
[[http://didgeridoo.une.edu.au/km/download.php?file=example1.tar.gz|Download Example1]]
One of the simplest analysis in quantitative genetics is the case where we have one trait of interest which is recorded only once, and where we have no additional random effects. Assuming a general pedigree (i.e. information from various types of relatives), we then fit a simple animal model, with additive genetic effects due to the individual to estimate the additive genetic and residual variance components and calculate the heritability of the trait from them.
==== Data and pedigree files ====
This is illustrated in Example 1 using some simulated data. The first few lines of the file ''tstmo1.d''
are shown in this box:
10025 1 2 1 220.0000
10026 1 2 1 212.0000
10027 1 2 1 221.0000
10028 1 2 1 207.0000
10029 1 2 1 218.0000
10030 1 2 1 201.0000
10031 1 2 1 214.0000
10032 1 2 1 229.0000
Column 1 contains the animal identity and columns 2 and 3 give the codes for its sire and dam, respectively. Note that all are numeric codes and that animal codes must be higher than both sire and dam codes. In this case,
we have a single fixed effect -- generation -- with 2 levels which is listed in column 4. Column 5 finally contains the trait value. As the data file has only one record per individual and as it contains the pedigree information for all individuals (i.e. there is no additional pedigree information on parents without records), we can use the same file as pedigree file.
==== Parameter file ====
The parameter file to analyse these data then looks like this (note that line numbers are NOT part of the parameter files and added here for ease of referencing only!):
COM Example 1 from DFREML : Simple univariate analysis
PED ../tstmo1.d
DAT ../tstmo1.d
animal
sire
dam
fixeff 9
record
end
ANAL UNI
MODEL
RAN animal NRM
FIX fixeff
TR record
END
VAR animal 1
60
VAR error 1
40
* Line 1 is a comment line (optional) -- what follows the identifier COM is written at the start of the summary files generated by WOMBAT.
* Line 2 identifies the pedigree file to be used. Here, it is held in the parent directory to the working directory (identified by ''../''). As no other option is given, this file is expected to have (at least) three columns, with columns 1, 2 and 3 containing the animal, sire and dam IDs, respectively.
Lines within a `block' are indented here -- this is for clarity (and recommended), but not actually required by the program.
* Line 3 specifies the data file to be used. It is the beginning of a `block' -- consisting of lines 4 to 8 -- which describes the content of the columns of the file, and is terminated by the ''END'' statement (line 9).\\ For columns to be `used' in the analysis, the name given is to be followed by a number (must be on integer) giving the maximum number of levels for this effect (this can be larger but should not be 10 or 100 times larger). For a random effect associated with a pedigree file (''animals'' in this case) this number be omitted (or given as ''0'') as the program will determine the correct number automatically.
* Line 10 specifies the type of analysis to be carried out: ''UNI'' for a standard, univariate analysis.
* Lines 11 to 15 specify the model of analysis. Again, this is a block of information (line 12 to 14) which is bracketed by the ''MODEL'' (line 11) and ''END'' (line 15) statements.
* Lines 16 to 19 finally give the variance components to be estimated and the starting values to be used in the iterative estimation scheme.
* There must be a line starting with ''VAR'' for each random effect fitted. This is followed by one or more lines with the starting values.
* In addition there must be ''VAR'' statement with the name ''error'' or ''residual'', and corresponding line(s) of starting values.
In this example, we call the parameter file ''wombat.par''.
==== Running WOMBAT ====
Running WOMBAT then requires typing in the program name, possibly with the appropriate path, together with any run time options desired on the `command line' Here we use
wombat -v
**N.B.** This requires a set-up where the operating system knows where to look for the executable file! Otherwise,
you may need to specify the complete path, e.g. ''/home/kmeyer/bin/wombat''. The run time option ''-v'' requests `verbose' output to the screen. Since our parameter file has the default name ''wombat.par'', we don't need to give its name.
==== Screen output ====
This generates the screen output as shown in the following:
Parameter file opened : "wombat.par"
**COM Example 1 from DFREML : Simple univariate analysis**
**PED ../tstmo1.d**
Pedigree file : "../tstmo1.d"
**DAT ../tstmo1.d**
Data file : "../tstmo1.d"
**animal**
Data - column : 1 "animal" -1
**sire**
Data - column : 2 "sire" -1
**dam**
Data - column : 3 "dam" -1
**fixeff 9**
Data - column : 4 "fixeff" 9
**record**
Data - column : 5 "record" -1
**end**
**ANAL UNI**
Analysis type = 1 "UNI"
No. of traits = 1
**MODEL**
**RAN animal NRM**
**FIX fixeff**
**TR record**
Trait no. 1 "record" Missing value -123456789
**END**
**VAR animal 1**
**60**
**VAR error 1**
**40**
End of parameter file reached : 19 lines read
Data file appers to be sorted correctly
No. of records = 282
No. of animal IDs in data = 282
No. of animal IDs in data = 282
Total no. of IDs in reduced ped file = 306
Pedigree file = ../tstmo1.d
Total no. of animal IDs found = 306
no. of NRM matrices 1 1
Pruned pedigree for random effect no. 1
Original no. of animals = 306
New no. of animals = 306
Inbreeding coefficients for random effect 1 computed
No. of inbred animals = 38
Average inbreeding coefficient = 2.2059 (in %)
... amongst inbred animals = 17.7632 (in %)
NRM Inverse no. : 1 "animal"
No. of elements = 906 log determinant = -195.46750491790496
No. of "subjects" in data 282
maxnr = 2 ncombi = 1
Parameters to be estimated & their starting values
1 COVS Z 1 1 40.000000
2 COVS A 1 1 60.000000
END data
Coefficient matrix fixed effects - rank = 2 no. of rows = 2
END "lsqfixed"
No. of equations in full rank submatrix = 308
Elements in MMM - after data 282
Elements in MMM - after RE 1 882
End of "GENMMD" 308 nofsub = 950
Elements in MMM - after RHS 1166
End of "SMBFCT": maxsub 318
maxlnz 1274
nnops 4599. 999999999999999.
END "eqnorders"
***log L 1 -760.345 Time 0
search 1
***AI iterate 0 log L = -760.344597 D = 0.0000 0.0000 Time 0
AI matrix : min. & max. eigenvalue 0.580946 6.50567 tau = 0.0000
***log L 2 -759.598 Time 0
Tau = 0.0000 Alpha = 1.0000 log L = -759.598
***AI iterate 1 log L = -759.597721 D = 0.74688 0.18539 Time 0
Estimates for "tau" = 0.0000 "alpha" = 1.0000 "norm" = 1.6662
1 CHOL Z 1 1 Old = 6.324555 New = 7.396408
2 CHOL A 1 1 Old = 7.745967 New = 6.298573
AI matrix : min. & max. eigenvalue 0.689232 6.18917 tau = 0.0000
***log L 3 -759.507 Time 0
Tau = 0.0000 Alpha = 1.0000 log L = -759.507
***AI iterate 2 log L = -759.506748 D = 0.90973E-01 0.52339E-01 Time 0
Estimates for "tau" = 0.0000 "alpha" = 1.0000 "norm" = 0.64462
1 CHOL Z 1 1 Old = 7.396408 New = 7.097490
2 CHOL A 1 1 Old = 6.298573 New = 6.713435
AI matrix : min. & max. eigenvalue 0.637744 6.31627 tau = 0.0000
***log L 4 -759.504 Time 0
Tau = 0.0000 Alpha = 1.0000 log L = -759.504
***AI iterate 3 log L = -759.503659 D = 0.30884E-02 0.11687E-01 Time 0
Estimates for "tau" = 0.0000 "alpha" = 1.0000 "norm" = 0.72789E-01
1 CHOL Z 1 1 Old = 7.097490 New = 7.148321
2 CHOL A 1 1 Old = 6.713435 New = 6.611619
AI matrix : min. & max. eigenvalue 0.652444 6.34003 tau = 0.0000
***log L 5 -759.503 Time 0
Tau = 0.0000 Alpha = 1.0000 log L = -759.503
***AI iterate 4 log L = -759.503454 D = 0.20554E-03 0.29391E-02 Time 0
Estimates for "tau" = 0.0000 "alpha" = 1.0000 "norm" = 0.19922E-01
1 CHOL Z 1 1 Old = 7.148321 New = 7.134470
2 CHOL A 1 1 Old = 6.611619 New = 6.636684
*** AIREML has converged
Last change in log L 2.05541655532215373E-4 5.0000000000000001E-4
1 animal 306 zero soln.s 0
SS residuals 9470.5877529900336
"WOMBAT" has finished !
==== Results files ====
=== Summary of pedigree information ===
The first output file generated is a summary of the pedigree information found. This file is named ''SumPedigree.out''.
======= Version 15-09-2009 ======================================= **KM** ====
Program WOMBAT : Summary of Pedigree Information
==============================================================================
Example 1 from DFREML : Simple univariate analysis
Analysis type : "UNI"
Data file : "../tstmo1.d"
Pedigree file : "../tstmo1.d"
Parameter file : "wombat.par"
* Lines 1 to 4 represent a header to describe the file. This includes the date (a.k.a version) that the program executable has been compiled (line 1).
* Line 6 contains the comment line given in the parameter file.
* Lines 8 to 11 identify the type of analysis to be carried out and the names of the data, pedigree and parameter files used.
* Lines 12 and 13 are blank.
No. of animal IDs in data file = = 282
No. of animal IDs in total = = 306
*****Pedigree Structure for random effect : 1 ****************************
Original no. of animals = 306
No. of animals after pruning = 306
... proportion (%) remaining = 100.0
No. of levels w/out records = 24
No. of levels with records = 282 100.0%
... 1 record(s) = 282 100.0%
No. of animals w/out offspring = 258 84.3%
No. of animals with offspring = 48 15.7%
... and records = 24 7.8%
No. of animals with unknown sire = 24
No. of animals with unknown dam = 24
No. of animals with both parents unknown = 24
No. of animals with records =
... and unknown sire = 0
... and unknown dam = 0
... and both parents unknown = 0
No. of sires = 12
... with progeny in the data = 12
... with records & progeny in data = 6
No. of dams = 36
... with progeny in the data = 36
... with records & progeny in data = 18
No. of animals with known/unpruned grand-parents
... with paternal grandsire = 144
... with paternal granddam = 144
... with maternal grandsire = 144
... with maternal granddam = 144
Inbreeding coefficients for random effect 1 computed
No. of inbred animals = 38
Average inbreeding coefficient = 2.2059 (in %)
... amongst inbred animals = 17.7632 (in %)
random effect no. = 1 NRM
no. of elements in NRM/GIN inverse 906
log determinant = -195.46750491790496
======== end of file ============================21-09-2009==========15:40====
* Lines 14 and 15 give the numbers of distinct animal codes (ID for identity) that WOMBAT has identified in the data and pedigree file, respectively.
* Lines 17 to 48 contain various counts on numbers of individuals and the numbers of ancestors identified together with some statistics on individuals which have both progeny and records in the data.
* Lines 50 to 56 report on the calculation of inbreeding coefficients and the characteristics of the inverse of the numerator relationship matrix that has been set up from the pedigree information supplied.
* Finally, line 57 marks the end of the file, together with a dat and time stamp.
=== Summary for data file and model specified ===
The second summary file generated at the beginning of a run is ''SumModel.out''. It gives a summary of the number of records, means and standard deviations per trait and the number of levels for each effect fitted that WOMBAT has found.
* It is a good idea to inspect this file closely `early on' to be sure that you and WOMBAT agree on the details of the analysis! WOMBAT provides the run-time option ''––setup'' to facilitate this.
======= Version 15-09-2009 ======================================= **KM** ====
Program WOMBAT : Summary of information from Set-up step
==============================================================================
Example 1 from DFREML : Simple univariate analysis
Analysis type : "UNI"
Data file : "../tstmo1.d"
Pedigree file : "../tstmo1.d"
Parameter file : "wombat.par"
No. of traits = 1
nrec mean sdev min. max.
1 "record" 282 229.787 13.4938 194.000 269.000
Fixed effects
1 "record" nlev
1 "fixeff" 2
Random effects nlev
1 "animal" 306 NRM
======== end of file ============================21-09-2009==========15:40====
=== Summary of variance component estimates ===
The main results file is named ''SumEstimates.out''.
======= Version 15-09-2009 ======================================= **KM** ====
Program WOMBAT : Estimates of covariance components
==============================================================================
Example 1 from DFREML : Simple univariate analysis
Analysis type : "UNI"
Data file : "../tstmo1.d"
Pedigree file : "../tstmo1.d"
Parameter file : "wombat.par"
No. of traits = 1 record
No. of records = 282 282
No. of parameters = 2
Maximum log L = -759.503
-1/2 AIC & AICC = -761.503 -761.525
-1/2 BIC = -765.138 "Penalty factor" = 2.817
Parameter estimates with approx. sampling erors
1 CHOL Z 1 1 7.13447 0.662468
2 CHOL A 1 1 6.63668 1.12245
Convergence criteria for last 3 iterates
Change in log likelihood = 0.090973 0.003088 0.000206
Change in parameter vector = 0.052339 0.011687 0.002939
Norm of gradient vector = 0.6446 0.0728 0.0199
Newton decrement = -0.2226 -0.0083 -0.0005
* Lines 1 to 11 again comprise the file header, together with information on the input files used.
* Lines 16 to 19 give the number of parameters estimated and the maximum of the REML log likelihood at the end of the run, together with corresponding information criteria calculated from it.
***** Estimates of residual covariances ************************************
Order of fit = 1
Covariance matrix
1 50.901
Matrix of correlations and variance ratios
1 0.5361
Covariances & correlations with approximate sampling errors
1 COVS Z 1 1 50.9007 9.45272 vrat 0.536 0.123
***** Estimates for RE 1 "animal" ***************************************
No. of levels = 306
Covariance structure = NRM
Order of fit = 1
Covariance matrix
1 44.046
Matrix of correlations and variance ratios
1 0.4639
Covariances & correlations with approximate sampling errors
2 COVS A 1 1 44.0456 14.8987 vrat 0.464 0.123
***** Estimates of phenotypic covariances ***********************************
Covariance matrix
1 94.946
Covariances & correlations with approximate sampling errors
3 COVS T 1 1 94.9462 10.0250
======== end of file ============================21-09-2009==========15:40====
==== Other output to look at ====