Re: coding R structure in a bivariate model

Dear Alistair,

re:

I am trying to fit a bivariate animal model where one trait (BW) has
repeated measures on 408 animals, (with phenotype recorded at 13 ages),
and the other trait (ELF) is measured only once. I've performed a
univariate analysis of BW (part 1 below) using a simple random
regression model (on age) for the additive effect, a pe effect, and a
diagonal error structure that has 13 age-specific residual variances.
The univariate analysis of ELF (part 2 below) is a straightforward
animal model partitioning variance into Va and Vr.

However, I'd like to run a bivariate model (part 3) and I'm having
trouble working out how to code the R structure. What I have below (part
3) works but fits a single Vr across all
age classes for BW. What I'd like to end up with the age specific Vr's
for BW (as in my univariate model), a Vr for ELF, and a COVr between BW
at each age and ELF, i.e. a structure like:

?
0 ?
0 0 ?
0 0 0 ?
0 0 0 0 ?
0 0 0 0 0 ?
0 0 0 0 0 0 ?
0 0 0 0 0 0 0 ?
0 0 0 0 0 0 0 0 ?
0 0 0 0 0 0 0 0 0 ?
0 0 0 0 0 0 0 0 0 0 ?
0 0 0 0 0 0 0 0 0 0 0 ?
0 0 0 0 0 0 0 0 0 0 0 0 ?
? ? ? ? ? ? ? ? ? ? ? ? ? ?

Where "?" is a (co)variance to be estimated
Below is the code I have so far, if anyone can offer any suggestions I'd
be very grateful.

Best,

Alastair

*** COMMENT1 : Since these are repeated measures, there needs to be
covariance
betwwen times as well as between traits.

####################################
!part 1 # random regression model for BW

# 12 age-specific records for each of 408 animals
#ageF is age as 13 level factor

BW ~ mu age !r !{animal animal.age !} ide(animal)

1 2 2 !STEP 0.001

408
ageF 0 DIAG !S2==1
1 1 1 1 1 1 1 1 1 1 1 1 1 #13 starting values

animal 2
animal 0 US
1.277
0.1132 0.1287
animal

ide(animal) 1
ide(animal) 0 ID 0.1

*** COMMENT 2. This analysis assumes the data is arranged sorted ageF
within animal.

*** COMMENT 3. This model has a common residual covariance between times
but
different variances at each time. Assuming the times represent growth
from near birth
and big changes in variance, the risk is that some of the 'DIAG' variances
will be negative
and the analysis fail. Since this has not happened, the common covariance
must have been small.

*** COMMENT 4 an alternative model to try is

BW ~ mu age !r !{animal animal.age !} pol(age,1).ide(animal)

1 2 2 !STEP 0.001

408
ageF 0 DIAG !S2==1
1 1 1 1 1 1 1 1 1 1 1 1 1 #13 starting values

animal 2
animal 0 US
1.277
0.1132 0.1287
animal

pol(age,1).ide(animal) 2
pol(age,1) 0 US !GP
0.1 0.01 0.01
ide(animal) 0 ID

*** COMMENT 5 From this model, you could test whether
a common variance for the DIAG values was plausible.
It is easiest to fit in the bivariate case.

###########################

!part 2 #univariate animal model for ELF

ELF ~ mu !r animal

*** Comment 6; for this to be correct on the same data file,
the ELF value must be either on a separate record, or on only
one data record per animal (i.e. not repeated on every animal record).

##############################

!part 3 #bivariate model

BW ELF ~ Trait at(Tr,1).age !r !{Trait.animal at(Tr,1).animal.age
!} at(Tr,1).ide(animal)

1 2 2

0
Trait 0 US !S2==1
1
0 1

Trait.animal 2
3 0 US
0.3
0.1 1.2
0.1 0.1 0.13
animal 0 AINV

at(Tr,1).ide(animal) 1
ide(animal) 0 ID 0.1

*** COMMENT 7.
Assuming, ELF is only present on one (say the first) of the BW records,

!part 3 #bivariate model

BW ELF ~ Trait at(Tr,1).age !r !{ Trait.animal at(Tr,1).animal.age !}
    !{ Trait.ide(animal) at(Tr,1).ide(animal).age !}

1 2 2

0
Trait 0 US !GPZF !S2==1
1
0 .01 # Fix this variance at say 1 hundredth of the real value
        # The other 99% dill be fitted in the G structure.

Trait.animal 2
3 0 US
0.3
0.1 1.2
0.1 0.1 0.13
animal 0 AINV

Trait.ide(animal) 2
3 0 US
0.3
0.1 1.2
0.1 0.1 0.13

ide(animal) 0 ID

*** COMMENT 8
We have artificially split the ELF residual into two components
and will get the error covariance from the larger component.

*** COMMENT 9
It is important to have the form of covariance at the genetic level
the same or simpler than at the residual level. Otherwise, the
variances may be wrongly partitioned.

I trust this helps. Let's know how well it works.

May Jesus Christ be gracious to you,

Arthur Gilmour, His servant .

Mixed model regression mapping for QTL detection in experimental crosses.
Computational Statistics and Data Analysis 51:3749-3764 now available at
http://dx.doi.org/10.1016/j.csda.2006.12.031

Profile: http://www.dpi.nsw.gov.au/reader/17263
Personal website: http://www.cargovale.com.au/

mailto:Arthur.Gilmour_at_dpi.nsw.gov.au, arthur_at_cargovale.com.au
Principal Research Scientist (Biometrics)
NSW Department of Primary Industries
Orange Agricultural Institute, Forest Rd, ORANGE, 2800, AUSTRALIA

fax: 02 6391 3899; 02 6391 3922 Australia +61
telephone work: 02 6391 3815; home: 02 6364 3288; mobile: 0438 251 426

ASREML 2 is now available from http://www.VSNi.co.uk/products/asreml
The ASReml discussion group is at ASREML-L_at_dpi.nsw.gov.au
To join it, mailto:arthur.gilmour_at_dpi.nsw.gov.au
Cookbook: http://uncronopio.org/ASReml

Proposed travel:
Leave early September
AAABG 24-26 September Armidale
    <><><><><><><><><><><><><><><><><><><><><><><><><>

This message is intended for the addressee named and may contain
confidential information. If you are not the intended recipient or
received it in error, please delete the message and notify sender. Views
expressed are those of the individual sender and are not necessarily the
views of their organisation.
Received on Tue Jun 15 2007 - 12:26:12 EST