Skip to content
Prev 858 / 20628 Next

user defined covariance structure in lme4 /nlme (Tapan Mehta)

Thanks Stephan. I appreciate your response. I agree with you about
convergence issue. I ran into this when trying on a slightly more complicated
pedigree data and had to use compound
symmetry matrix.

 I have tried using corSymm but it does not give me correct random
effects estimates as compared to what I get through SAS. I am trying to
estimate the genetic and error variance and based on a simulated data
(through SAS) with genetic variance 7 and error variance 2 whereas the
results I obtain from the code below gives 7.58 and 1.38. 

migGdata<-c(0.0,0.5,0.5)
cs1Symm<-corSymm(migGdata,form=~1 | FAM_ID)
cs2Symm<-Initialize(cs1Symm,data=migTrio)
test<-lme(y ~ x1 + x2, random=~1|FAM_ID,correlation=cs2Symm,data=migTrio)

Let G = k*C, where G is the variance of the random effects, C is a
symmetric positive definite matrix, and k a parameter. We would like to
supply lme /lmer  with C (user defined covariance structure) and let it
estimate k. I am also trying to use the pdSymm method to implement this
but I am running into this error:

 > mt
     [,1] [,2] [,3]
[1,]  1.0  0.0  0.5
[2,]  0.0  1.0  0.5
[3,]  0.5  0.5  1.0
pd2 <- pdSymm(mt,form=~1|FAM_ID,name=c("1","2","3"))
testRes<-lme(Y~X1+X2,random=pd2,data=migTrio)

Error in getGroups.data.frame(dataMix, groups) : 
Invalid formula for groups

Please let me know if any of you can suggest what I need to fix to solve this problem.

Thanks,

Tapan


----- Original Message ----
From: Stephan Moratti <moratti at med.ucm.es>
To: r-sig-mixed-models at r-project.org
Sent: Wednesday, April 16, 2008 3:05:20 AM
Subject: [R-sig-ME] user defined covariance structure in lme4 /nlme (Tapan Mehta)

Hi Tapan,

The following line using lme from the nlme package should work:

testRes <- lme(Y ~X1+X2,random=~1|FAM_ID,correlation=corSymm(form = ~1|FAM_ID),data=fam_dat)

testRes <- lme(Y ~X1+X2,random=~1|FAM_ID,correlation=corSymm(),data=fam_dat) should be equal.

I hope this helps. Sometimes I had the problem that the REML did not converge using the corSymm(), whereas using a compound symmetry structure worked.

Best,

Stephan