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

matthew s06mw3 at abdn.ac.uk
Wed Dec 4 19:15:48 CET 2013


Dear Fiona,

You can check out the course notes to the MCMCglmm package - there are
examples in there that simulate data.  Also, try the "grfx" and "drfx"
functions in the "nadiv" package.  Below is some code using your example
that does what I think you want.  Of course, there are probably other
(more elegant) options out there that can be accessed by a little
googling...

Sincerely,
Matthew

######## R CODE   ######
library(MCMCglmm)
library(nadiv)
# make a random dataset
data<-data.frame(group=factor(rep(seq(50),100)))
#define the mean for each trait
mus <- rnorm(3, 100, 20)
#define the desired variance-covariance matrix for the "groups"
vcov.in <- matrix(c(10,  5, 3,
                      5, 10, 8,
              3,  8, 10), nrow = 3, byrow = TRUE)
#define a residual variance-covariance matrix (here, residuals are
un-correlated among groups)
Rvcov <- diag(c(10, 20, 5))
#Create the random effects for the group variable
fx.out <- drfx(G = vcov.in, fac = "group", dataf = data)
#create the residuals
r.out <- grfx(n = 5000, G = Rvcov, incidence = NULL)
#create a trait value for each of the three traits
data$x <- mus[1] + fx.out$fx[, 1] + r.out[, 1]
data$y <- mus[2] + fx.out$fx[, 2] + r.out[, 2]
data$z <- mus[3] + fx.out$fx[, 3] + r.out[, 3]
# 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)

#now compare the covariance matrix of generated random effects with the
estimated one from MCMCglmm
cov(fx.out$fx)
varcov



--
....................................................
Dr. Matthew E. Wolak
School of Biological Sciences
Zoology Building
University of Aberdeen
Tillydrone Avenue
Aberdeen AB24 2TZ
office phone: +44 (0)1224 273255




The University of Aberdeen is a charity registered in Scotland, No SC013683.



More information about the R-sig-mixed-models mailing list