[R-sig-ME] simulating datasets for mixed model analysis

[Fiona Ingleby]

Dear list,

I am trying to simulate a dataset with a specific covariance matrix for analysis with a mixed model. To illustrate what I mean, I've used MCMCglmm (although a similar model could be done using lmer) to make a reproducible example:

# make a random dataset
data<-data.frame(x=rnorm(5000,0.1,0.1),y=rnorm(5000,0.5,0.2),z=rnorm(5000,0.4,0.1),group=factor(rep(seq(50),100)))
# analyse with MCMCglmm as if x, y and z are three response variables and group is a grouping factor
prior<-list(R=list(V=diag(3)/3,nu=0.05),G=list(G1=list(V=diag(3)/3,nu=0.05)))
model<-MCMCglmm(cbind(x,y,z)~trait-1,random=~us(trait):group,
              rcov=~us(trait):units,prior=prior,data=data,family=rep("gaussian",3),
              nitt=20000,burnin=1000,thin=25)
# extract the group-level variance/covariance matrix:
varcov<-matrix(c(mean(model$VCV[,1]),mean(model$VCV[,2]),mean(model$VCV[,3]),mean(model$VCV[,4]),mean(model$VCV[,5]),mean(model$VCV[,6]),mean(model$VCV[,7]),mean(model$VCV[,8]),mean(model$VCV[,9])),3)

What I would like to do is start with a given matrix, and simulate a dataset around this matrix - so, something similar to the line of code above used to create the random dataset, but designed for a known var/cov structure.

I'd be really grateful if anyone could offer any help or point me in the right direction at all.

Thanks in advance,

Fiona