Print

Bivariate models

Summary

Using bivariate models allows the estimation of parameters relating to each trait alone (i.e. VA, h2 etc), but also yields estimates of covariance components between traits. These include the (additive) genetic covariance COVA, which is often rescaled to give the genetic correlation rG. However, correlations among traits can also arise through other random effects (e.g. maternal effects) and these sources can be explicitly modelled in a bivariate analysis as well.

More details on bivariate models including sample data sets and code are available in the tutorial available here.

Table of contents



Definitions

Sample Code

ASReml

The code for a multivariate model is similar to the univariate case, but a few extra lines are required to specify the model of the (co)variance structures we want to fit. These extra lines are the principal source of confusion for new ASReml users but are necessary since the program can actually fit a wide variety of strutures. The simplest - an unstructured covariance matrix - is often appropriate and is used in the code below.

Given a data set where two traits (BWT and TARSUS) are measured in a set of individuals across a number of years, we can partition variance for each trit into components arising from addiive effects, year effect, and residual (unexplained) effects. We can also apply the same partitioning to the covariance between traits. Note that "Trait" is the multivariate equivalent of mu (i.e. fits a mean for each trait)

ASReml bivariate animal model
  ANIMAL  !P
  BWT
  TARSUS
  SEX     !A  #sex as a factor
  YEAR    !A  # year as a factor

pedigreedata.ped      !skip 1   
phenotypicdata.dat    !skip 1    

BWT TARSUS ~ Trait Trait.SEX !r Trait.ANIMAL  Trait.YEAR

1 2 2
0
Trait 0 US     !S2==1   
0.1
0.01 0.1

Trait.ANIMAL 2
Trait 0   US
0.1               #starting values for VA1
0.01 0.1          #                    COVA  VA2
ANIMAL  0 !AINV

Trait.YEAR 2
Trait 0   US
0.1
0.01 0.1
YEAR


If we wanted to test the significance of the genetic covariance here (against a null hypothesis that COVA=0) we could compare the log-likelihood of this model to one in which we modified the covariance matix associated with ANIMAL such that COVA is constrained to equal zero (there is one additional parameter and hence one less DF in the unconstrained model).
Trait.ANIMAL 2
Trait 0  US !GPZP  #forces VA1 and VA2 to be positive (P) but the COVA to be zero.
0.1
0 0.1
ANIMAL  0 !AINV


ASReml-R

To run a multivariate analysis in ASReml-R you have to use cbind to bind your response variables together. To fit an intercept for each trait, you have to use 'trait' as the intercept and to fit the fixed effects for both variables, interact the effect with 'trait'. In a bivariate model for each random effect you will have three outputs -the variance component for each response variable and the covariance between the two.

The following code fits a bivariate model of birth wt and tarsus length, with intercept and sex as fixed effects and additive genetic effects and year as random effects.
#create modela, see numbered comments below for description of each line 
modela<-asreml(cbind(BWT,TARSUS)~trait +trait:SEX #1 
  ,random=~us(trait,init=c(1,0.1,1)):ped(ANIMAL)  #2 
          +us(trait,init=c(1,0.1,1)):year
  ,rcov= ~ units:us(trait,init=c(0.1,0.1,0.1))    #3 
  ,data=gryphon                                   #4 
  ,ginverse=list(ANIMAL=ainv)                     #5 
  ,maxiter=20)                                    #6 
#1: Bivariate model of BWT & TARSUS, fitting the mean of each and sex as fixed effects
#2: Random G structure defined as an unstructured covariance matrix (us). Initial 
starting values of va1=1 cov12=0.1 va2=1 
#3: Residual (R) structure defined as an unstructured covar matrix, initial starting 
values of va1, cov12, va2. 
#4: Use data file gryphon 
#5: connect the individual in the data file to the pedigree 
#6: maximum iterations =20

Testing significance of a covariance

To test the significance of a covariance, fix the value of the covariance to zero and then compare the models with and without the covariance using log-likelihood ratio tests. In ASReml-R, you can do this by specifying the covariance matrix as a diagonal matrix (i.e. diag instead of us). To test the significance of the genetic covariance in the above model, use the following code, note the the number of starting values has also decreased.
modela<-asreml(cbind(BWT,TARSUS)~trait +trait:SEX #1 
  ,random=~us(trait,diag=c(1,1)):ped(ANIMAL)      #2 
          +us(trait,init=c(1,0.1,1)):year
  ,rcov= ~ units:us(trait,init=c(0.1,0.1,0.1))    #3 
  ,data=gryphon                                   #4 
  ,ginverse=list(ANIMAL=ainv)                     #5 
  ,maxiter=20)                                    #6


MCMCglmm

Fitting a bivariate model requires a little more coding. The critical changes relative to the univariate models are (1) we need to provide priors in the form of matrices (2) we use the cbind() to combine arrays of response variable observations, (3) we use the reserved variable 'trait' to specify fixed and random effects, (4) we must now explicitly set up a model specification section to handle the residual covariance, and (5) we need to supply MCMCglmm with the type of distribution to fit for each response variable. Additionally, it will generally be necessary to run multivariate models for longer than univariate models in order to ensure adequate mixing. Here we have coded 10 times the default run parameters.

# observed phenotypic variance for priors
phen.var<-matrix(c(var(Data$BWT,na.rm=TRUE),0,0,
                        var(Data$TARSUS,na.rm=TRUE)),2,2)

# an example prior specification for a bivariate
# model with a single random effect
prior2.1<-list(G=list(G1=list(V=phen.var/2,n=2)),
                           R=list(V=phen.var/2,n=2))

# a simple bivariate model with only a genetic effect
model2.1<-MCMCglmm(cbind(BWT,TARSUS)~trait-1,
                                 random=~us(trait):animal,
                                 rcov=~us(trait):units,
                                 family=c("gaussian","gaussian"),
                                 pedigree=Ped,data=Data,
                                 nitt=130000,thin=100,burnin=30000,
                                 prior=prior2.1,verbose=FALSE)

# validity checks
plot(model2.1$VCV[,"animal.trait.TARSUS.TARSUS"])

autocorr(model2.1$VCV)[,,"animal.trait.TARSUS.TARSUS"][3,4]

posterior.mode(model2.1$VCV)


Heritabilities and genetic correlations, as well as their statistical support are easily obtained:

# h2=Va/Vp
heritability.BWT2.1<-model2.1$VCV[,"animal.trait.BWT.BWT"]/
                         (model2.1$VCV[,"animal.trait.BWT.BWT"]+
                             model2.1$VCV[,"units.trait.BWT.BWT"])

posterior.mode(heritability.BWT2.1)

heritability.TARSUS2.1<-model2.1$VCV[,"animal.trait.TARSUS.TARSUS"]/
                         (model2.1$VCV[,"animal.trait.TARSUS.TARSUS"]+
                            model2.1$VCV[,"units.trait.TARSUS.TARSUS"])

posterior.mode(heritability.TARSUS2.1)

# CORRg=COVa/sqrt(Va1*Va2)
genetic.correlation2.1<-model2.1$VCV[,"animal.trait.BWT.TARSUS"]/
                          sqrt(model2.1$VCV[,"animal.trait.BWT.BWT"]*
                              model2.1$VCV[,"animal.trait.TARSUS.TARSUS"])

posterior.mode(genetic.correlation2.1)


Fitting additional random effects is similar to the univariate case, only we need to continue specifying them in the same manner as the genetic effect in the multivariate context, i.e. using the reserved variable 'trait'.

# an example prior specification for a model
# with three random effects
prior2.2<-list(G=list(G1=list(V=phen.var/4,n=2),
                      G2=list(V=phen.var/4,n=2),
                      G3=list(V=phen.var/4,n=2)),
                      R=list(V=phen.var/4,n=2))

# an example bivariate model with three random effects
model2.2<-MCMCglmm(cbind(BWT,TARSUS)~trait-1+trait:SEX,
            random=~us(trait):animal+us(trait):BYEAR+us(trait):MOTHER,
            rcov=~us(trait):units,
            family=c("gaussian","gaussian"),
            pedigree=Ped,data=Data,
            nitt=130000,thin=100,burnin=30000,
            prior=prior2.2,verbose=FALSE)

# very basic validity check
autocorr(model2.2$VCV)[,,"animal.trait.TARSUS.TARSUS"][3,4]

# example post-processing for the most 
# interesting/interpretable parameters
posterior.mode(model2.2$VCV)

genetic.correlation2.2<-model2.2$VCV[,"animal.trait.BWT.TARSUS"]/
                sqrt(model2.2$VCV[,"animal.trait.BWT.BWT"]*
                   model2.2$VCV[,"animal.trait.TARSUS.TARSUS"])

maternal.correlation2.2<-model2.2$VCV[,"MOTHER.trait.BWT.TARSUS"]/
                sqrt(model2.2$VCV[,"MOTHER.trait.BWT.BWT"]*
                   model2.2$VCV[,"MOTHER.trait.TARSUS.TARSUS"])

posterior.mode(genetic.correlation2.2)

posterior.mode(maternal.correlation2.2)

HPDinterval(genetic.correlation2.2,0.95)

HPDinterval(maternal.correlation2.2,0.95)


WOMBAT

Running multivariate analyses in WOMBAT requires only slightly more complex code than a univariate analysis. To run a bivariate analysis, write the following command file using your text editor. The main differences are in the parts describing the data structure, the model structure and the starting values.
# A single optional comment line of up to 74 characters
COMMENT WOMBAT analysis of Gryphon birth weight and tarsus length

# The type of analysis to be performed 
# In this case a MUltiVariate analysis with 2 traits
ANALYSIS MUV 2

# Name the pedigree file, assuming it is in the same folder as the parameter file
PEDS gryphon.ped

# Name the data file, assuming it is in the same folder as the parameter file
# GRP indicates that we will be using the 'compact' data specification (see manual)
DATA gryphon_bi.dat GRP

  TRNOS 1 2
# Each variable to be fitted in the model needs to be followed by
# the maximum number of levels
  TRAITNO
  ANIMAL 900
  MOTHER 400
  BYEAR 35
  SEX 2
  NAMES BWT TARSUS
END

# Model specification
# Type of effect (FIXed, COVariate, RANdom, TRAIT) and variable names
# NRM indicates that a pedigree is available for ANIMAL
MODEL
  RAN ANIMAL NRM
  TRAIT BWT 1
  TRAIT TARSUS 2
END MODEL

# Specify starting values as upper triangle of variance-covariance matrix
# and the number of rows and columns in the matrix
# It doesn't matter if the values are next to each other,
# below each other, or something in between
VAR  residual 2
1.0
0.1
1.0
VAR  ANIMAL 2
1.0 0.1 1.0

Note that the starting values supplied here are arbitrary. If the model is difficult to fit then it can be because the starting values are too far from the best estimates. One way around this is to run single trait models first to get good starting values for the variances (but you still have to “guess” starting values for the covariances).

UPDATE: Significance testing of genetic correlations

To formally test the significance of COVA we would have to compare the log-likelihood for this model to one in which we specify that COVA = 0. Unfortunately, at this moment it is not possible to constrain individual covariance components in WOMBAT. This is currently one of its main limitations. It is however possible to constrain all the off-diagonal elements of the model to zero, which in a bivariate model comes down to the same thing. In this case we only have to specify starting values for the diagonal elements (i.e. the variances) for the ANIMAL term. We do this by specifying the starting values as follows:
VAR  residual 2
1.0
0.1
1.0
VAR  ANIMAL 2 DIAG
1.0
1.0
VAR  BYEAR 2
1.0 0.1
    1.0
VAR  MOTHER 2
1.0 0.1 1.0



Menu

Main Menu

Recently visited pages