[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