[R] animal models and lme

[Jarrod Hadfield]

-- 
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.98747062,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.11508781,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.27469486,-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.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,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.125,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.125,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.125,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.125,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.25,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.125),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.25,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