Post processing

Introduction

ASReml includes a procedure to calculate certain functions of variance components either as a final stage of an analysis or as a post-analysis procedure. These functions enable the calculation of heritabilities and correlations from simple variance components and when US, CORUH and XFA structures are used in the model fitting. A simple example is shown in the code box. The instructions to perform the required operations are listed after the VPREDICT !DEFINE line and terminated by a blank line. ASReml holds the instructions in a .pin until the end of the job when it retrieves the relevant information from the .asr and .vvp files and performs the specified operations. The results are reported in the .pvc file.

 VPREDICT !DEFINE
 F PhenVar Sire + Residual
 F GenVar Sire * 4
 H herit GenVar PhenVar

Syntax

Instructions to calculate functions are headed by a line VPREDICT !DEFINE. This line and the following instructions can occur anywhere in the .as file but the logical place is at the end of the file. The instructions are processed after the job (part/cycle) has been completed. ASReml recognises a blank line (or end of file) as termination of the functional instructions.

Functions of the variance components are specified by lines of the form
letter label coefficients

letter (either C, D, F, H, R, S, V or X) must occur in column 1

C is a copy function to facilitate rearranging components,

D is a delete function,

F forms linear combinations of variance components,

H is for forming heritabilities, the ratio of two components,

R is for forming the correlation from a covariance component,

S is a square root function,

V is for converting components related to a CORUH or an XFA structure into components related to a US structure,

X is a multiply function,

label names the result,

coefficients is the list of arguments/coefficients for the linear function.

When ASReml reads back the variance parameters from the .asr file, each component is assigned a name and a sequential number. The name is the model term, possibly qualified by the component, and may relate to a set of parameters. Names can be abbreviated to a unique substring. For example the set of parameters \textcolor{blue} us(Trait).nrm(Tag);us(Trait) may be referenced as Tag;us. Individual parameters within a set are specified using subscripts in square braces as in Tag;us[1].

The original implementation was based entirely on the numbers but it will generally be better to use the names, since the order model terms are reported cannot always be predicted.

Linear combinations of components

First ASReml extracts the variance components from the .asr file and their variance matrix from the .vvp file.

The F, S, V and X functions create new components which are appended to the list. For example, the F function appends component k + c'v and forms Cov(c'v, v) and Var(c'v) where v is the vector of existing variance components, c is the vector of coefficients for the linear function and k is an optional offset which is usually omitted but would be 1 to represent the residual variance in a probit analysis and 3.289 to represent the residual variance in a logit analysis. The general form of the directive is
F label a + b * c_b+ c + d + m*k
where a, b, c and d are the numbers or names of existing components v_a, v_b, v_c and v_d and c_b is a multiplier for v_b
m is a number bigger than the current length of v to flag the special case of adding the offset k. m*k is usually omitted but k would be 1 to represent the residual variance in a probit analysis and 3.289 to represent the residual variance in a logit analysis.

When using the component numbers, the form a:b can be used to reference blocks of components as in
F label a:b * k + c:d

 VPREDICT !DEFINE
 F phenvar 1 + 2 # Sire + Res
 F genvar 1 * 4  # Sire * 4
 H herit 4 3     # genvar phenvar

Assuming that the instructions in the ASReml code box corresponds to a simple sire model ( ~ mu !r Sire) so that variance component 1 is the Sire variance and variance component 2 is the Residual variance, then

F phenvar 1 + 2	or
F phenvar Sire + Residual	creates a third component called phenvar
	which is the sum of the variance components,
	that is, the phenotypic variance,
F genvar 1 * 4	or
F genvar Sire * 4	creates a fourth component called genvar
	which is the sire variance component multiplied by 4,
	that is, the genotypic variance, and
H herit 4 3	or
H herit genvar phenvar	calculates the heritability.

Copy and Delete components

C copy i[:j]=k[:l] copies the component k (and following) to i (and following).

D discard i resets the number of components held (back) to i-1.

Convert CORUH and XFA to US

S label i:j when i:j are assumed positive variance parameters, inserts components which are the SQRT of components i:j,

V label i:j where i:j spans an XFA variance structure, inserts the US matrix based on the XFA parameters,

V label i:j where i:j spans an CORUH variance structure, inserts the US matrix based on the CORUH parameters,

X label i*k inserts a component being the product of components i and k

X label i:j*k inserts j-i+1 components being the products of components i:j and k

X label i:j*k:l inserts a set of j-i+1 components being the pairwise products of components i:j and k:l

The C, D, S, V and X functions are new in ASReml 4. The multiply option (X) allows a correlation in a CORUV structure to be converted to a covariance. The SQRT option allows conversion of CORGH to US, provided the dimension is moderate (say <10).

The basic operations used in deriving the variances are:
Cov(v_A,v_Av_B) = v_ACov(v_A,v_B) + v_BVar(v_A)
Var(v_Av_B) = v_A²Var(v_B) + 2v_Av_BCov(v_A,v_B) + v_B²Var(v_A)
Cov(v_C,v_Av_B)=v_ACov(v_B,v_C)+v_BCov(v_A,v_C)
Var(v_A+v_B)=Var(v_A)+2Cov(v_A,v_B)+Var(v_B)
Cov(v_A,v_A+v_B) = Var(v_A)+Cov(v_A,v_B)
Cov(v_C,v_A+v_B)= Cov(v_A,v_C)+Cov(v_B,v_C)
Var(√ v_A) = Var(v_A)/(4v_A)
Cov(v_C,√ v_A) = Cov(v_A,v_C)/(2√ v_A)
Cov(v_A,√ v_A) = Var(v_A)/(2√ v_A)

Heritability

Heritabilities are requested by function lines beginning with an H. The specific form of the directive is
H label n d
This calculates σ²_n/σ²_d and se[σ²_n/σ²_d] where n and d are the names of the components or integers pointing to components v_n and v_d that are to be used as the numerator and denominator respectively in the heritability calculation.

Var(σ²_n/σ²_d) =(σ²_n/ σ²_d)² ( Var(σ²_n)/ (σ²_nσ²_n) + Var(σ²_d)/ (σ²_dσ²_d) - 2 Cov(σ²_n, σ²_d)/ (σ²_nσ²_d))

 VPREDICT !DEFINE
 F phenvar 1 + 2 # Sire + Res
 F genvar 1 * 4  # Sire * 4
 H herit 4 3     # genvar phenvar

In the example

H herit 4 3	or
H herit genvar phenvar	calculates the heritability by calculating
	component 4 (from second line) / component 3 (from first line),
	that is, `genetic variance / phenotypic variance`.

Correlation

Correlations are requested by lines beginning with an R. The specific form of the directive is
R label a ab b
This calculates the correlation r = σ_ab/√(σ²_aσ²_b) and the associated standard error. a, b and ab are integers indicating the position of the components to be used. Alternatively,
R label a:n
calculates the correlation r = σ_ab/√(σ²_aσ²_b) for all correlations in the lower triangular row-wise matrix represented by components a to n and the associated standard errors. }

 VPREDICT !DEFINE
 F phenvar 1:3 + 4:6  #F phenvar Sire + Residual
 R phencorr 7 8 9     #R phencorr phenvar
 R gencorr 4:6        #R gencorr Sire

Var(r)=r²[Var(σ²_a) / ( 4 σ²_aσ²_a) +Var(σ²_b) / (4 σ²_bσ²_b) +Var(σ_ab} / σ²_ab
+2Cov(σ²_a,σ²_b)/ (4 σ²_aσ²_b) -2Cov(σ²_a,σ_ab ) / (2 σ²_aσ_ab ) -2Cov(σ_ab,σ²_b) / (2 σ_abσ²_b)]

In the example

R phencorr 7 8 9	or
R phencorr phenvar	calculates the phenotypic covariance by calculating
	`v₈ / √( v₇v₉)}` where components 7, 8 and 9 are created
	with the first line of the .pin file, and
R gencorr 4:6	or
R gencorr Sire	calculates the genotypic covariance by calculating
	`v₅ / √(v₄v₆)` where components 4, 5 and 6
	are variance components from the analysis.

A more detailed example

The following example for a {em bivariate sire model} is a little more complicated. The job file bsiremod.as contains

 ...
 coop.fmt

 ywt fat ~ Trait Trait.(age c(brr) sex sex.age) !r us(Trait).sire  !f Tr.grp
 residual units.us(Trait)

 VPREDICT !DEFINE
 F phenvar units.us(Trait);us(Trait) + us(Trait).sire;us(Trait) # 1:3 + 4:6
 F addvar sire * 4                                              # 4:6 * 4
 H heritA addvar[1] phenvar[1]                                  # 10 7
 H heritB addvar[3] phenvar[3]                                  # 12 9
 R phencorr phenvar                                             # 7 8 9
 R gencorr addvar                                               # 4:6

The relevant lines of the .asr file are

 Model_Term                             Sigma         Sigma   Sigma/SE   % C
 units.us(Trait)               1356 effects
 Trait                   US_V  1  1   19.7352       19.7352      18.02   0 P
 Trait                   US_C  2  1   2.58252       2.58252      10.37   0 P
 Trait                   US_V  2  2   1.70532       1.70532      18.04   0 P
 us(Trait).sire                 184 effects
 Trait                   US_V  1  1   1.78978       1.78978       2.09   0 P
 Trait                   US_C  2  1  0.573268E-01  0.573268E-01   0.43   0 P
 Trait                   US_V  2  2  0.582773E-01  0.582773E-01   1.54   0 P

Numbering the parameters reported in bsiremod.asr (and bsiremod.vvp)
v₁ error variance for ywt
v₂ error covariance for ywt and fat
v₃ error variance for fat
v₄ sire variance component for ywt
v₅ sire covariance for ywt and fat
v₆ sire variance for fat
then

F phenvar units.us(Trait);us(Trait) + us(Trait).sire;us(Trait)	or
F phenvar units + sire	or
F phenvar 1:3 + 4:6
	creates new components `v₇= v₁+ v₄,`
	v₈= v₂+ v₅ and `v₉= v₃+ v₆`,
F addvar sire * 4	or
F addvar 4:6 * 4	creates new components v₁₀ `=` 4 v₄
	v₁₁ `=` 4 v₅ and v₁₂ `=` 4v₆,
H heritA addvar[1] phenvar[1]	or
H heritA 10 7	forms v₁₀ / v₇ to give the heritability for ywt,
H heritB addvar[3] phenvar[3]	or
H heritB 12 9	forms v₁₂ / v₉ to give the heritability for fat,
R phencorr phenvar	or
R phencorr 7 8 9	forms v₈`/√(v₇v₉)`,
	that is, the phenotypic correlation between ywt and fat,
R gencorr addvar	or
R gencorr 4:6	forms v₅/√(v₄v₆),
	that is, the genetic correlation between ywt and fat.

The resulting .pvc file contains:

           - - - Results from analysis of ywt fat - - -

 units.us(Trait)               1356 effects
   1 units.us(Trait);us(Trait)               V  1  1      19.7352         1.09518
   2 units.us(Trait);us(Trait)               C  2  1      2.58252        0.249038
   3 units.us(Trait);us(Trait)               V  2  2      1.70532        0.945299E-01
 us(Trait).sire                 184 effects
   4 us(Trait).sire;us(Trait)                V  1  1      1.78978        0.856354
   5 us(Trait).sire;us(Trait)                C  2  1     0.573268E-01    0.133318
   6 us(Trait).sire;us(Trait)                V  2  2     0.582773E-01    0.378424E-01
   7 phenvar  1                21.525        1.3639
   8 phenvar  2                2.6398       0.27639
   9 phenvar  3                1.7636       0.99638E-01
  10 addvar  4                 7.1591        3.4274
  11 addvar  5                0.22931       0.53416
  12 addvar  6                0.23311       0.15106
     heritA       = addvar    10/phenvar    7=          0.3326    0.1476
     heritB       = addvar    12/phenvar    9=          0.1322    0.0836
     phenco  2  1 = phenv  8/SQR[phenv  7*phenv  9]=    0.4285    0.0347
     gencor  2  1 = addva 11/SQR[addva 10*addva 12]=    0.1775    0.3835
 Notice: The parameter estimates are followed by
          their approximate standard errors.

The first 6 lines are copied from the .asr file.

VPREDICT: PIN file processing

There are four forms of the VPREDICT directive.

If the .pin file exists and has the same name as the jobname (including any suffix appended by using !RENAME), just specify the VPREDICT directive.

If the .pin file exists but has a different name to the jobname, specify the VPREDICT directive with the .pin file name as its argument.

If the .pin file does not exist or must be reformed, a name argument for the file is optional but the !DEFINE qualifier should be set. Then the lines of the .pin file should follow on the next lines, terminated by a blank line.

An alternative to using VPREDICT is process the contents of the .pin file by running ASReml with the -P command line option specifying the .pin file as the input file.

Return to index