[R] animal models and lme

[Simon Blomberg]

Hi,

You should look at Pinheiro and Bates (2000) Mixed-effects models in S and S-Plus. It describes how to format the correlation matrix to pass to functions lme and gls. Basically, the correlation matrix has to be one of the corStruct classes, probably corSymm for your example. So in the call to lme (or gls if you really have no random effects), use something like: correlation=corSymm(AGRM[lower.tri(AGRM)], fixed=TRUE). Anyway, it's all there in the book.

HTH,

Simon.

Simon Blomberg, PhD
Depression & Anxiety Consumer Research Unit
Centre for Mental Health Research
Australian National University
http://www.anu.edu.au/cmhr/
Simon.Blomberg at anu.edu.au  +61 (2) 6125 3379


> -----Original Message-----
> From: Jarrod Hadfield [mailto:jarrod.hadfield at imperial.ac.uk]
> Sent: Wednesday, 23 July 2003 5:52 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] animal models and lme
> 
> 
> 
> -- 
> Hi R users,
> 
> I'm trying to fit an animal model using lme and cannot fit the 
> corelation matrix which describes the additive genetic correlation 
> between individuals.
> 
> I have created a test data set (pedigree) where the first column 
> specifies the id of each individual, and the second column contains 
> simulated trait values. AGRM is the corresponding additive genetic 
> relationship matrix whose elements are equal to the genetic 
> correlation between individuals.
> 
> The test data is simulated to have a heritibility of 0.8 (i.e. 80% of 
> the variation can be explained by genetic effects)
> 
> My initial attempt has been:
> 
> CM<-pdSymm(AGRM, ~1|animal)
> lme(trait~1, data=pedigree, CM)
> 
> but which ever combination I use I always get "invalid 
> formula for groups".
> 
> Does any one know how to code this sort of model?
> 
> Thanks for your help,
> 
> Jarrod Hadfield.
> 
> copy and paste from here to get the test data (Sorry the correlation 
> matrix is so big)
> 
> ##############################################################
> ################################
> 
> animal<-c("Newbird1","Newbird2","a","b","c","d","e","f","g","h
> ","i","j","k","l","m","n","o","p","q","r","s","t","u","v","w",
> "x","y","z","aa","bb","cc","dd","ee","ff","gg","hh","ii","jj",
> "kk","ll","mm")
> 
> trait<-c(-1.10988343,-0.95061690,-1.05596407,-0.20075671,-0.98
747062,1.72132706,1.33973268,0.30333354,0.34969961> ,0.70681745,-0.15210172,-1.05918719,-0.32665059,0.55423824,-0.
> 81324407,-0.61359451,,0.65936390,-0.92295239,0.87285686,1.1150
> 8781,1.30904325,0.83437615,,0.72948628,,0.63043694,-0.40078431
> ,,1.20814060,,0.69233357,-1.51568342,-0.16235203,-2.26109886,1
.21242313,0.04010594,0.06988673,-0.70644384,-> 0.69870290,-0.02043294,-0.57957532,-0.78726879,0.20574982,0.27
> 469486,-0.25899545)
> 
> pedigree<-data.frame(cbind(animal, trait))
> 
> AGRM<-as.matrix(cbind(c(1,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.5,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0),c(0,1,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.5,0.5,0.5),c(0,0,1,0,0,0,0,0,0,0,0.5,0.5,0.
> 5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.25,0.25,0.25,0,0,0
> ,0,0,0,0,0,0),c(0,0,0,1,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,0.25,0.25,0.25,0,0,0,0,0,0,0,0,0),c(0,0,0,
0,1,0,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0,0,0,0,0,0,> 0,0,0,0,0,0.5,0.5,0.5,0.25,0.25,0.25,0.25,0.25,0.25,0,0,0),c(0
,0,0,0,0,1,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0,0,0,0> ,0,0,0,0,0,0,0,0,0,0,0.75,0.75,0.75,0.25,0.25,0.25,0,0,0),c(0,
0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,> 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),c(0,0,0,0,0,0,0,1,0,0,0
> ,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
> ,0,0,0),c(0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,
> 0.5,0.5,0,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.25,0.25,0.25),c(0,0,0,
0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5!
>  , 
> 0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0.25,0.25),c(0,0,0.5,0.5,
> 0,0,0,0,0,0,1,0.5,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0.25,0.
> 25,0.25,0.25,0.25,0,0,0,0,0,0,0,0,0),c(0,0,0.5,0.5,0,0,0,0,0,0
> ,0.5,1,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0.25,0.25,0.25,0.2
> 5,0.25,0,0,0,0,0,0,0,0,0),c(0,0,0.5,0.5,0,0,0,0,0,0,0.5,0.5,1,
> 0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0.25,0.25,0.25,0.25,0.25,0,0,
> 0,0,0,0,0,0,0),c(0,0,0.5,0.5,0,0,0,0,0,0,0.5,0.5,0.5,1,0,0,0,0
> ,0,0,0,0,0,0,0,0,0.25,0.25,0.25,0.5,0.5,0.5,0,0,0,0,0,0,0,0,0)
> ,c(0,0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,1,0.5,0.5,0.5,0,0,0,0,0,0,0
> ,0,0,0,0,0.25,0.25,0.25,0.75,0.75,0.75,0.25,0.25,0.25,0,0,0),c
> (0,0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,0.5,1,0.5,0.5,0,0,0,0,0,0,0,0
> ,0,0,0,0.25,0.25,0.25,0.5,0.5,0.5,0.25,0.25,0.25,0,0,0),c(0,0,
0,0,0.5,0.5,0,0,0,0,0,0,0,0,0.5,0.5,1,0.5,0,0,0,0,> 0,0,0,0,0,0,0,0.25,0.25,0.25,0.5,0.5,0.5,0.25,0.25,0.25,0,0,0)
> ,c(0,0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,0.5,0.5,0.5,1,0,0,0,0,0,0,0
> ,0,0,0,0,0.25,0.25,0.25,0.5,0.5,0.5,0.5,0.5,0.5,0,0,0),c(0,0,0
> ,0,0,0,0.!
>  5 
> ,0.5,0,0,0,0,0,0,0,0,0,0,1,0.5,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0
> ,0,0,0,0,0,0,0),c(0,0,0,0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0.5,
> 1,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),c(0,0,0,0,0,0
,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0.5,0.5,1,0.5,0,0,0,0> ,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),c(0,0,0,0,0,0,0.5,0.5,0,0,0,0,
> 0,0,0,0,0,0,0.5,0.5,0.5,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
> 0),c(0,0,0,0,0,0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,1,0.5,0.5
> ,0.5,0,0,0,0,0,0,0,0,0,0.25,0.25,0.25,0.25,0.25,0.25),c(0,0,0,
0,0,0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.5,1,0.> 5,0.5,0,0,0,0,0,0,0,0,0,0.25,0.25,0.25,0.25,0.25,0.25),c(0,0,0
,0,0,0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5> ,1,0.5,0,0,0,0,0,0,0,0,0,0.25,0.25,0.25,0.5,0.5,0.5),c(0,0,0,0
,0,0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,0> .5,1,0,0,0,0,0,0,0,0,0,0.25,0.25,0.25,0.25,0.25,0.25),c(0.5,0,
0.5,0,0,0,0,0,0,0,0.25,0.25,0.25,0.25,0,0,0,0,0,0,> 0,0,0,0,0,0,1,0.5,0.5,0.125,0.125,0.125,0,0,0,0,0,0,0,0,0),c(0
.5,0,0.5,0,0,0,0,0,0,0,0.25,0.25,0.25,0.25,0,0,0,0> ,0,0,0,!
>  0 
> ,0,0,0,0,0.5,1,0.5,0.125,0.125,0.125,0,0,0,0,0,0,0,0,0),c(0.5,
0,0.5,0,0,0,0,0,0,0,0.25,0.25,0.25,0.25,0,0,0,0,0,> 0,0,0,0,0,0,0,0.5,0.5,1,0.125,0.125,0.125,0,0,0,0,0,0,0,0,0),c
> (0,0,0.25,0.25,0.5,0,0,0,0,0,0.25,0.25,0.25,0.5,0.25,0.25,0.25
> ,0.25,0,0,0,0,0,0,0,0,0.125,0.125,0.125,1,0.5,0.5,0.125,0.125,
> 0.125,0.125,0.125,0.125,0,0,0),c(0,0,0.25,0.25,0.5,0,0,0,0,0,0
> .25,0.25,0.25,0.5,0.25,0.25,0.25,0.25,0,0,0,0,0,0,0,0,0.125,0.
> 125,0.125,0.5,1,0.5,0.125,0.125,0.125,0.125,0.125,0.125,0,0,0)
> ,c(0,0,0.25,0.25,0.5,0,0,0,0,0,0.25,0.25,0.25,0.5,0.25,0.25,0.
> 25,0.25,0,0,0,0,0,0,0,0,0.125,0.125,0.125,0.5,0.5,1,0.125,0.12
> 5,0.125,0.125,0.125,0.125,0,0,0),c(0,0,0,0,0.25,0.75,0,0,0,0,0
> ,0,0,0,0.75,0.5,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0.125,0.125,0.12
> 5,1.25,0.75,0.75,0.25,0.25,0.25,0,0,0),c(0,0,0,0,0.25,0.75,0,0
> ,0,0,0,0,0,0,0.75,0.5,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0.125,0.12
> 5,0.125,0.75,1.25,0.75,0.25,0.25,0.25,0,0,0),c(0,0,0,0,0.25,0.
> 75,0,0,0,0,0,0,0,0,0.75,0.5,0.5,0.5,0,0,0,0,0,0,0,0,0,0,0,0.12
> 5,0.125,0!
>  . 
> 125,0.75,0.75,1.25,0.25,0.25,0.25,0,0,0),c(0,0,0,0,0.25,0.25,0
> ,0,0.5,0,0,0,0,0,0.25,0.25,0.25,0.5,0,0,0,0,0.25,0.25,0.25,0.2
> 5,0,0,0,0.125,0.125,0.125,0.25,0.25,0.25,1,0.5,0.5,0.125,0.125
> ,0.125),c(0,0,0,0,0.25,0.25,0,0,0.5,0,0,0,0,0,0.25,0.25,0.25,0
> .5,0,0,0,0,0.25,0.25,0.25,0.25,0,0,0,0.125,0.125,0.125,0.25,0.
> 25,0.25,0.5,1,0.5,0.125,0.125,0.125),c(0,0,0,0,0.25,0.25,0,0,0
> .5,0,0,0,0,0,0.25,0.25,0.25,0.5,0,0,0,0,0.25,0.25,0.25,0.25,0,
> 0,0,0.125,0.125,0.125,0.25,0.25,0.25,0.5,0.5,1,0.125,0.125,0.1
> 25),c(0,0.5,0,0,0,0,0,0,0.25,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25
> ,0.25,0.5,0.25,0,0,0,0,0,0,0,0,0,0.125,0.125,0.125,1,0.5,0.5),
> c(0,0.5,0,0,0,0,0,0,0.25,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0.2
> 5,0.5,0.25,0,0,0,0,0,0,0,0,0,0.125,0.125,0.125,0.5,1,0.5),c(0,
0.5,0,0,0,0,0,0,0.25,0.25,0,0,0,0,0,0,0,0,0,0,0,0,> 0.25,0.25,0.5,0.25,0,0,0,0,0,0,0,0,0,0.125,0.125,0.125,0.5,0.5,1)))
> colnames(AGRM)<-animal
> rownames(AGRM)<-animal
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>