Multivariate animal genetics data - Sheep

Introduction

The analysis of incomplete or unbalanced multivariate data often presents computational difficulties. These difficulties are exacerbated by either the number of random effects in the linear mixed model, the number of traits, the complexity of the variance models being fitted to the random effects or the size of the problem. In this section we illustrate two approaches to the analysis of a complex set of incomplete multivariate data.

Much of the difficulty in conducting such analyses in ASReml centres on obtaining good starting values. Derivative based algorithms such as the AI algorithm can be unreliable when fitting complex variance structures unless good starting values are available. Poor starting values may result in divergence of the algorithm or slow convergence. A particular problem with fitting unstructured variance models is keeping the estimated variance matrix positive definite. These are not simple issues and in the following we present a pragmatic approach to them.

The data are taken from a large genetic study on Coopworth lambs. A total of 5 traits, namely weaning weight ( wwt), yearling weight ( ywt), greasy fleece weight ( gfw), fibre diameter ( fdm) and ultrasound fat depth at the C site ( fat) were measured on 7043 lambs. The lambs were the progeny of 92 sires and 3561 dams, produced from 4871 litters over 49 flock-year combinations. Not all traits were measured on each group. No pedigree data was available for either sires or dams.

The aim of the analysis is to estimate heritability (h2) of each trait and to estimate the genetic correlations between the five traits. We will present two approaches, a half-sib analysis and an analysis based on the use of an animal model, which directly defines the genetic covariance between the progeny and sires and dams.

The data fields included factors defining sire, dam and lamb ( tag), covariates such as age, the age of the lamb at a set time, brr the birth rearing rank (1 = born single raised single, 2 = born twin raised single, 3 = born twin raised twin and 4 = other), sex (M, F) and grp a factor indicating the flock-year combination.

Half-sib analysis

In the half-sib analysis we include terms for the random effects of sires, dams and litters. In univariate analyses the variance component for sires is denoted by σ2s= 0.25σ2A where σ2A is the additive genetic variance, the variance component for dams is denoted by σ2d= 0.25σ2A+ σ2m where σ2m is the maternal variance component and the variance component for litters is denoted by σ2l and represents variation attributable to the particular mating. For a multivariate analysis these variance components for sires, dams and litters are, in theory replaced by unstructured matrices, one for each term. Additionally we assume the residuals for each trait may be correlated. Thus for this example we would like to fit a total of 4 unstructured variance models. For such a situation, it is sensible to commence the modelling process with a series of univariate analyses. These give starting values for the diagonals of the variance matrices, but also indicate what variance components are estimable. The ASReml job for the univariate analyses is 
 Multivariate Sire  Dam model
   tag
   sire   92 !I
   dam  3561 !I
   grp    49
   sex
   brr     4
   litter 4871
   age       wwt   !M0 ywt   !M0    # !M0 recodes zeros as missing values
   gfw   !M0 fdm   !M0 fat   !M0
 coop.fmt
 wwt ~ mu age brr sex age.sex !r sire dam lit age.grp sex.grp !f grp
Table 1. REML estimates of a subset of the variance parameters for each trait for the genetic example, expressed as a ratio to their asymptotic s.e.
term wwt ywt gfw fdm fat
sire 3.68 3.57 3.95 1.92 1.92
dam 6.25 4.93 2.78 0.37 0.05
litter 8.79 0.99 2.23 1.91 0.00
age.grp 2.29 1.39 0.31 1.15 1.74
sex.grp 2.90 3.43 3.70 - 1.83

Wald tests of the fixed effects for each trait for the genetic example,
term wwt ywt gfw fdm fat
age 331.3 67.1 52.4 2.6 7.5
brr 554.6 73.4 14.9 0.3 13.9
sex 196.1 123.3 0.2 2.9 0.6
age.sex 10.3 1.7 1.9 - 5.0

Tables 1 and 2 present the summary of these analyses. Fibre diameter was measured on only 2 female lambs and so interactions with sex were not fitted. The dam variance component was quite small for both fibre diameter and fat. The REML estimate of the variance component associated with litters was effectively zero for fat.

Thus in the multivariate analysis we consider fitting the following models to the sire, dam and litter effects,

Var(us) = Σs ⊗ I92
Var(ud) = Σd ⊗ I3561
Var(ul) = Σl ⊗ I4891
where Σs (5 by 5), Σd (3 by 3) and Σl (4 by 4) are positive definite symmetric matrices corresponding to the between traits variance matrices for sires, dams and litters respectively. The variance matrix for dams does not involve fibre diameter and fat depth, while the variance matrix for litters does not involve fat depth. The effects in each of the above vectors are ordered levels within traits. Lastly we assume that the residual variance matrix is given by
Σe ⊗ I7043

Table 3 presents the sequence variance models fitted to each of the four random terms sire, dam, litter and error in the ASReml job 
 !WORK 1 !CONTINUE  !RENAME !ARG 1 2 3
 Multivariate Sire  Dam model
   tag
   sire   92 !I
   dam  3561 !I
   grp    49
   sex
   brr     4
   litter 4871
   age       wwt   !m0 ywt   !m0    # !M0 identifies missing values
   gfw   !m0 fdm   !m0 fat   !m0
 coop.fmt   !DOPATH $1  !MAXIT 20
 !SUBSET DamTr Trait 1 2 3 0 0
 !SUBSET LitTr Trait 1 2 3 4 0

  wwt ywt gfw fdm fat ~ Trait Tr.age Tr.brr Tr.sex Tr.age.sex,
 !PATH 1 // !r  xfa1(Tr).sire xfa1(DamTr).dam xfa1(LitTr).litter,

 !PATH 2 // !r  us(Tr).sire  xfa1(DamTr).dam us(LitTr).litter,

 !PATH 0 // at(Trait,1,2,4,5).age.grp    at(Trait,1,2,3,5).sex.grp  !f Tr.grp
 residual units.us(Trait)
Table 3. Variance models fitted for each part of the ASReml job in the analysis of the genetic example
term matrix !PATH 1 !PATH 2
sire Σs FA1 US
dam Σd FA1 XFA1
litter Σl FA1 US
error Σe US US

In !PATH 1, the error variance model is taken to be unstructured using the starting values are from the sample covariance matrix of the data. For incomplete data the matrix so obtained may not, in general be positive definite. The command to run !PATH 1 is 
 asreml -n mt 1
The Loglikelihood from this run is -20000-1426.18. The results from this analysis can be automatically used by ASReml for the next part, if the next part is run immediately, or if the .rsv is copied prior to running the next part. That is, we add the !PATH 2 coding to the job, copy mt1.rsv to mt2.rsv so that when we run !PATH 2 it starts from where !PATH 1 finished, and run the job using 
 asreml -cnrw1 mt 2
The Loglikelihood from this run is -20000-1419.99.

ASReml has 3 forms of fitting the factor analytic model. The fa() form uses the correlation scale and really only works for a single FA factor. The xfa() form used here is generally preferred. Note the 4 fitted matrices for the sire, dam and litter components: 
 Covariance/Variance/Correlation Matrix US Residual
   9.485      0.5689      0.2375      0.1494      0.2172
   7.361       17.65      0.4317      0.1596      0.4618
  0.2766      0.6859      0.1430      0.3499      0.1665
  0.8655       1.261      0.2489       3.538      0.4268E-01
  0.8300       2.407      0.7815E-01  0.9963E-01   1.540
 Covariance/Variance/Correlation Matrix XFA xfa1(Tr).sire
  0.5808      0.6878      0.4472E-01 -0.6127E-01  0.4116      0.7061
  0.6515       1.545      0.6170E-01 -0.8452E-01  0.5678      0.9741
  0.4204E-02  0.9459E-02  0.1522E-01 -0.5496E-02  0.3692E-01  0.6334E-01
 -0.1976E-01 -0.4446E-01 -0.2869E-03  0.1791     -0.5058E-01 -0.8677E-01
  0.5808E-01  0.1307      0.8433E-03 -0.3963E-02  0.3428E-01  0.5830
  0.5381       1.211      0.7813E-02 -0.3672E-01  0.1079      1.0000
 Covariance/Variance/Correlation Matrix XFA xfa1(DamTr).dam
   2.168      0.9377      0.8010      0.9377
   2.296       2.765      0.8542      1.0000
  0.1736      0.2091      0.2167E-01  0.8542
   1.381       1.663      0.1257      1.0000
 Covariance/Variance/Correlation Matrix XFA xfa1(LitTr).litter
   3.508      0.5418     -0.2023      0.2412E-01 -0.9972
   1.417       1.950     -0.1102      0.1314E-01 -0.5433
 -0.4601E-01 -0.1869E-01  0.1474E-01 -0.4908E-02  0.2029
  0.4528E-01  0.1840E-01 -0.5972E-03   1.004     -0.2419E-01
  -1.868     -0.7588      0.2463E-01 -0.2424E-01  1.0000
The hign correlations among dam effects for the traits indicate that it may not be possible to fit an unstructured variance matrix for this term. Indeed, ASReml cannot form a positive definite unstructured dam matrix. The dam components are substantially larger than the sire components indicating a large mjaternal environment effect on wwt, ywt and gfw. Migrating the tag and litter terms to unstructured only increased the LogL by 6.19 but added 7 variance parameters, so the fit is not significantly improved.

Note that by running these models in order, and using !CONTINUE, ASReml uses the parameter estimates from the previous model as starting values for the next model. A portion of the final output file is 
  11 LogL=-1419.99     S2= 1.00000      18085 df
  12 LogL=-1419.99     S2= 1.00000      18085 df

 Covariance/Variance/Correlation Matrix US Residual
   9.463      0.5690      0.2356      0.1639      0.2183
   7.329       17.54      0.4242      0.2498      0.4639
  0.2728      0.6688      0.1417      0.3995      0.1680
  0.9625       1.997      0.2870       3.643      0.4881E-01
  0.8333       2.411      0.7851E-01  0.1156       1.541
 Covariance/Variance/Correlation Matrix US us(Tr).sire
  0.5941      0.7042      0.2965      0.2016      0.2698
  0.6740       1.542      0.1333E-01 -0.1250      0.5727
  0.2798E-01  0.2027E-02  0.1499E-01  0.1069     -0.4610E-02
  0.6185E-01 -0.6180E-01  0.5211E-02  0.1585     -0.6350
  0.3783E-01  0.1293     -0.1027E-03 -0.4598E-01  0.3308E-01
 Covariance/Variance/Correlation Matrix XFA at(Tr,1).dam
   2.154      1.0000      0.8016      1.0000
   2.210       2.268      0.8016      1.0000
  0.1606      0.1648      0.1864E-01  0.8016
   1.468       1.506      0.1094      1.0000
 Covariance/Variance/Correlation Matrix US at(Tr,1).lit
   3.555      0.5098     -0.1174     -0.4092E-01
   1.542       2.574      0.2020     -0.5249
 -0.3059E-01  0.4479E-01  0.1910E-01 -0.3201
 -0.7311E-01 -0.7981     -0.4192E-01  0.8981
Note the XFA matrix have an extra (fourth) row/column. This last row/column pertains to the 'factor' and relates to the other rows via the factor 'loadings'.

In the .res file is reported an eigen analysis of these four variance structures. 
 Eigen Analysis of US matrix for Residual
 Eigen values    22.456     5.213     3.394     1.160     0.103
   Percentage    69.468    16.126    10.501     3.588     0.318
      1          0.4970   -0.8664    0.0147    0.0469    0.0027
      2          0.8509    0.4764   -0.1321   -0.1745   -0.0327
      3          0.0335    0.0230    0.0585   -0.0047    0.9974
      4          0.1170    0.0878    0.9842    0.0771   -0.0633
      5          0.1187    0.1194   -0.1013    0.9805    0.0039

 Eigen Analysis of US matrix for us(Tr).sire
 Eigen values     1.902     0.304     0.114     0.012     0.010
   Percentage    81.190    12.978     4.855     0.533     0.444
      1          0.4580    0.7479    0.4687   -0.1059    0.0065
      2          0.8859   -0.3645   -0.2770    0.0269   -0.0694
      3          0.0077    0.0794    0.0846    0.9506   -0.2880
      4         -0.0170    0.5256   -0.8015    0.1074    0.2635
      5          0.0710   -0.1588    0.2323    0.2701    0.9180

 Eigen Analysis of XFA matrix for xfa1(DamTr).dam
 Eigen values     5.434     0.007    -0.000    -0.000
   Percentage    99.878     0.122    -0.000    -0.000
      1          0.6296    0.0296    0.4766    0.4766
      2          0.6461    0.0304   -0.7590   -0.7590
      3          0.0470   -0.9989   -0.0000   -0.0000
      4          0.4290    0.0202    0.4437    0.4437

 Eigen Analysis of US matrix for at(Tr,1).lit
 Eigen values     4.758     1.808     0.464     0.016
   Percentage    67.523    25.662     6.583     0.233
      1          0.7840    0.5799    0.2207    0.0169
      2          0.6048   -0.6336   -0.4823   -0.0173
      3          0.0019   -0.0377    0.0161    0.9992
      4         -0.1399    0.5107   -0.8476    0.0333

Estimating genetic parameters

The REML estimates of all the variance matrices except for the dam components are positive definite. Heritabilities for each trait can be calculated using the VPREDICT/ .pin file facility of ASReml. The heritability is given by
h2 = σ2A/σ2P
where σ2P is the phenotypic variance given by
σ2P= σ2s+ σ2d+ σ2l+ σ2e
recalling that
σ2s = 0.25 σ2A
σ2d = 0.25 σ2A + σ2m

In the half-sib analysis we only use the estimate of additive genetic variance from the sire variance component. The ASReml .pin file is presented below along with the output from the following command 
 asreml -p mt2


 V DamVC  dam;xfa1   # Convert XFA parameters to variance components
 F AddVar sire;us * 4   # Convert sire variance to genetic variance
 F Phen units[1:6] + sire[1:6] + DamVC + litter[1:6]
 F Phen units[7:10] + sire[7:10] + litter[7:10]
 F Phen units[11:15] + sire[11:15]
 R PhenCor Phen
 R GenCor AddVar
 H HeritWWT AddVar[1] Phen[1]
 H HeritYWT AddVar[3] Phen[3]
 H HeritGFW AddVar[6] Phen[6]
 H HeritFDM AddVar[10] Phen[10]
 H HeritFAT AddVar[15] Phen[15]


 ASReml 4.1 [01 Apr 2014] Multivariate Sire Dam model
          mtsire02.pvc created 21 Oct 2014 14:23:22.023

          - - - Results from analysis of wwt ywt gfw fdm fat - - -

   1 Trait_1.age.grp                                     0.135181E-02    0.669213E-03
   2 Trait_2.age.grp                                     0.101813E-02    0.821073E-03
   3 Trait_4.age.grp                                     0.176601E-02    0.156284E-02
   4 Trait_5.age.grp                                     0.209288E-03    0.124576E-03
   5 Trait_1.sex.grp                                     0.920264        0.318430
   6 Trait_2.sex.grp                                      15.3787         4.39391
   7 Trait_3.sex.grp                                     0.279368        0.753013E-01
   8 Trait_5.sex.grp                                      1.44050        0.800278
 units.us(Trait)              35200 effects
   9 units.us(Trait);us(Trait)               V  1  1      9.46254        0.284160
  10 units.us(Trait);us(Trait)               C  2  1      7.32955        0.356496
  11 units.us(Trait);us(Trait)               V  2  2      17.5378        0.646913
  12 units.us(Trait);us(Trait)               C  3  1     0.272792        0.325139E-01
  13 units.us(Trait);us(Trait)               C  3  2     0.668806        0.476358E-01
  14 units.us(Trait);us(Trait)               V  3  3     0.141722        0.597479E-02
  15 units.us(Trait);us(Trait)               C  4  1     0.962553        0.333063
  16 units.us(Trait);us(Trait)               C  4  2      1.99706        0.547140
  17 units.us(Trait);us(Trait)               C  4  3     0.287033        0.565026E-01
  18 units.us(Trait);us(Trait)               V  4  4      3.64285        0.404761
  19 units.us(Trait);us(Trait)               C  5  1     0.833307        0.968962E-01
  20 units.us(Trait);us(Trait)               C  5  2      2.41105        0.124089
  21 units.us(Trait);us(Trait)               C  5  3     0.785061E-01    0.107543E-01
  22 units.us(Trait);us(Trait)               C  5  4     0.115625        0.947746E-01
  23 units.us(Trait);us(Trait)               V  5  5      1.54054        0.469248E-01
 us(Tr).sire                    460 effects
  24 us(Tr).sire;us(Tr)                      V  1  1     0.594134        0.161449
  25 us(Tr).sire;us(Tr)                      C  2  1     0.674027        0.211958
  26 us(Tr).sire;us(Tr)                      V  2  2      1.54206        0.396416
  27 us(Tr).sire;us(Tr)                      C  3  1     0.279826E-01    0.182893E-01
  28 us(Tr).sire;us(Tr)                      C  3  2     0.202626E-02    0.289466E-01
  29 us(Tr).sire;us(Tr)                      V  3  3     0.149924E-01    0.373875E-02
  30 us(Tr).sire;us(Tr)                      C  4  1     0.620564E-01    0.110815
  31 us(Tr).sire;us(Tr)                      C  4  2    -0.615623E-01   -0.380000
  32 us(Tr).sire;us(Tr)                      C  4  3     0.526314E-02    0.194931E-01
  33 us(Tr).sire;us(Tr)                      V  4  4     0.158546        0.861663E-01
  34 us(Tr).sire;us(Tr)                      C  5  1     0.378333E-01    0.406810E-01
  35 us(Tr).sire;us(Tr)                      C  5  2     0.129334        0.663251E-01
  36 us(Tr).sire;us(Tr)                      C  5  3    -0.103158E-03   -0.200000E-01
  37 us(Tr).sire;us(Tr)                      C  5  4    -0.459500E-01    -1.51000
  38 us(Tr).sire;us(Tr)                      V  5  5     0.330789E-01    0.160577E-01
 xfa1(DamTr).dam              14244 effects
  39 xfa1(DamTr).dam;xfa1(DamTr)             V  0  1     0.559130E-02    0.559130
  40 xfa1(DamTr).dam;xfa1(DamTr)             V  0  2      0.00000         0.00000
  41 xfa1(DamTr).dam;xfa1(DamTr)             V  0  3     0.666005E-02    0.528575E-02
  42 xfa1(DamTr).dam;xfa1(DamTr)             L  1  1      1.46524        0.180671
  43 xfa1(DamTr).dam;xfa1(DamTr)             L  1  2      1.50712        0.207022
  44 xfa1(DamTr).dam;xfa1(DamTr)             L  1  3     0.109589        0.216579E-01
 us(LitTr).litter             19484 effects
  45 us(LitTr).litter;us(LitTr)              V  1  1      3.55525        0.415333
  46 us(LitTr).litter;us(LitTr)              C  2  1      1.54283        0.464708
  47 us(LitTr).litter;us(LitTr)              V  2  2      2.56955        0.805502
  48 us(LitTr).litter;us(LitTr)              C  3  1    -0.305637E-01   -0.720000
  49 us(LitTr).litter;us(LitTr)              C  3  2     0.444459E-01    0.608848E-01
  50 us(LitTr).litter;us(LitTr)              V  3  3     0.190719E-01    0.781635E-02
  51 us(LitTr).litter;us(LitTr)              C  4  1    -0.731186E-01   -0.220000
  52 us(LitTr).litter;us(LitTr)              C  4  2    -0.798163        -1.56000
  53 us(LitTr).litter;us(LitTr)              C  4  3    -0.419291E-01   -0.760000
  54 us(LitTr).litter;us(LitTr)              V  4  4     0.898030        0.392153
  55 DamVC                     2.1525       0.33431
  56 DamVC                     2.2083       0.36993
  57 DamVC                     2.2714       0.62438
  58 DamVC                    0.16057       0.32569E-01
  59 DamVC                    0.16516       0.46275E-01
  60 DamVC                    0.18670E-01   0.58989E-02
  61 AddVar 24                 2.3765       0.64602
  62 AddVar 25                 2.6961       0.84894
  63 AddVar 26                 6.1682        1.5843
  64 AddVar 27                0.11193       0.73304E-01
  65 AddVar 28                0.81050E-02   0.11051
  66 AddVar 29                0.59970E-01   0.14970E-01
  67 AddVar 30                0.24823       0.44511
  68 AddVar 31               -0.24625       0.64456
  69 AddVar 32                0.21053E-01   0.76801E-01
  70 AddVar 33                0.63418       0.34523
  71 AddVar 34                0.15133       0.16341
  72 AddVar 35                0.51734       0.26553
  73 AddVar 36               -0.41263E-03   0.23627E-01
  74 AddVar 37               -0.18380       0.12146
  75 AddVar 38                0.13232       0.64267E-01
  76 Phen  9                   15.764       0.31272
  77 Phen 10                   11.755       0.37473
  78 Phen 11                   23.921       0.63109
  79 Phen 12                  0.43079       0.33010E-01
  80 Phen 13                  0.88044       0.44371E-01
  81 Phen 14                  0.19446       0.54925E-02
  82 Phen 15                  0.95149       0.29837
  83 Phen 16                   1.1373       0.37654
  84 Phen 17                  0.25037       0.37270E-01
  85 Phen 18                   4.6994       0.22536
  86 Phen 19                  0.87114       0.10440
  87 Phen 20                   2.5404       0.13825
  88 Phen 21                  0.78403E-01   0.12041E-01
  89 Phen 22                  0.69675E-01   0.97927E-01
  90 Phen 23                   1.5736       0.48632E-01
     PhenCo  2  1 = Phen  77/SQR[Phen  76*Phen  78]=    0.6053    0.0110
     PhenCo  3  1 = Phen  79/SQR[Phen  76*Phen  81]=    0.2460    0.0174
     PhenCo  3  2 = Phen  80/SQR[Phen  78*Phen  81]=    0.4082    0.0175
     PhenCo  4  1 = Phen  82/SQR[Phen  76*Phen  85]=    0.1105    0.0343
     PhenCo  4  2 = Phen  83/SQR[Phen  78*Phen  85]=    0.1073    0.0352
     PhenCo  4  3 = Phen  84/SQR[Phen  81*Phen  85]=    0.2619    0.0369
     PhenCo  5  1 = Phen  86/SQR[Phen  76*Phen  90]=    0.1749    0.0203
     PhenCo  5  2 = Phen  87/SQR[Phen  78*Phen  90]=    0.4141    0.0180
     PhenCo  5  3 = Phen  88/SQR[Phen  81*Phen  90]=    0.1417    0.0214
     PhenCo  5  4 = Phen  89/SQR[Phen  85*Phen  90]=    0.0256    0.0360
     GenCor  2  1 = AddVa 62/SQR[AddVa 61*AddVa 63]=    0.7042    0.1026
     GenCor  3  1 = AddVa 64/SQR[AddVa 61*AddVa 66]=    0.2965    0.1720
     GenCor  3  2 = AddVa 65/SQR[AddVa 63*AddVa 66]=    0.0133    0.1811
     GenCor  4  1 = AddVa 67/SQR[AddVa 61*AddVa 70]=    0.2022    0.3515
     GenCor  4  2 = AddVa 68/SQR[AddVa 63*AddVa 70]=   -0.1245    0.3249
     GenCor  4  3 = AddVa 69/SQR[AddVa 66*AddVa 70]=    0.1080    0.3871
     GenCor  5  1 = AddVa 71/SQR[AddVa 61*AddVa 75]=    0.2699    0.2725
     GenCor  5  2 = AddVa 72/SQR[AddVa 63*AddVa 75]=    0.5726    0.2023
     GenCor  5  3 = AddVa 73/SQR[AddVa 66*AddVa 75]=   -0.0046    0.2653
     GenCor  5  4 = AddVa 74/SQR[AddVa 70*AddVa 75]=   -0.6345    0.3773
     HeritWWT     = AddVar 2  61/Phen  9   76=          0.1508    0.0396
     HeritYWT     = AddVar 2  63/Phen 11   78=          0.2579    0.0624
     HeritGFW     = AddVar 2  66/Phen 14   81=          0.3084    0.0716
     HeritFDM     = AddVar 3  70/Phen 18   85=          0.1349    0.0717
     HeritFAT     = AddVar 3  75/Phen 23   90=          0.0841    0.0402
 Notice: The parameter estimates are followed by
          their approximate standard errors.
Note that the maternal genetic correlations have values much larger than 1. A better approach would be to use the animal model where we conclude there is no maternal effect for YWT and GFW.
  • Back

    Return to index