Solved – Multivariate model in lme() with independent random effect, similar to MCMCglmm

by

I would like to specify a multivariate model with lme with a random effect for group which is independent across variables. I found this post,
which explains that the model specified as:

fit.multilevel <- lme( y ~ var -1,  dd,                     random = ~ var -1| school,                        correlation =  corSymm( form = ~ v |school/id),                      weights = varIdent(form = ~ 1 | v))

is equivalent to (using MCMCglmm):

fit.w.null <- MCMCglmm( cbind(mathach, ses) ~ trait -1 ,         random = ~us(trait):school,            rcov = ~us(trait):units,         data = hs,         family = c("gaussian","gaussian"))

However the model I would really like has independent random effects for school across the variables:

fit.w.null <- MCMCglmm( cbind(mathach, ses) ~ trait -1 ,         random = ~idh(trait):school,            rcov = ~us(trait):units,         data = hs,         family = c("gaussian","gaussian"))

(the random effect is specified with idh()). I cannot find a way to specify this in lme.
Here is some code to generate example data and run the analysis:

gr=rep(c("A","B","C"),each=100) grp.mean=rep(rnorm(3,5,1),each=100) var1=rnorm(300,rep(rnorm(3,5,1),each=100),1) var2=0.8*var1+3*rnorm(300,rep(rnorm(3,5,1),each=100),2) dat.mcmc=data.frame(var1=var1,var2=var2,group=gr) dat.lme=data.frame(y=c(var1,var2),var=rep(c("1","2"),each=length(gr)),              group=rep(gr,2),id=rep(1:length(gr),2)) dat.lme$v=as.integer(dat.lme$var)

Currently this is the best I can do:

fit.multilevel <- lme( y ~ var -1,  dat.lme,                        random = ~ var -1| group,                           correlation =  corSymm( form = ~ v |group/id),                         weights = varIdent(form = ~ 1 | v))

This is (I believe) equivalent to:

fit.w.null <- MCMCglmm( cbind(var1, var2) ~ trait -1 ,                         random = ~us(trait):group,                            rcov = ~us(trait):units,                         data = dat.mcmc,                         family = c("gaussian","gaussian"))

However, I would like an lme command equivalent to:

fit.w.null <- MCMCglmm( cbind(var1, var2) ~ trait -1 ,                         random = ~idh(trait):group,                            rcov = ~us(trait):units,                         data = hs,                         family = c("gaussian","gaussian"))

Any help you can give me will be much appreciated.

Contents hide
1 Best Answer
2 Similar Posts:

Best Answer

I believe the answer is to use ?pdClasses as explained here.

fit.multilevel <- lme( y ~ var -1,  dat.lme,                    random = list(group= pdDiag(~ var-1)),                        correlation =  corSymm( form = ~ v |group/id),                     weights = varIdent(form = ~ 1 | v))

Similar Posts:

Rate this post

Leave a Comment