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

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")) ``

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))  ``