Hello, I would like to simulate nested data, where my mixed effects model fitted to real data has the form:
y ~ time + (1 | site/subject) I then take the hyper-parameters from this model to simulate fake data, using this function: create_fake <- function(J,K,L,HP,t){ # J : number of sites # K : number of subjects / site # L : number of years# HP: hyperparameters from fit, y ~ time + (1 | site/subject)# t: fractional effectiveness of treatment time <- rep(seq(0,2,length=L), J*K) subject <- rep(1:(J*K), each=L) site <- sample(rep (1:J, K)) site1 <- factor(site[subject]) treatment <- sample(rep (0:1, J*K/2)) treatment1 <- treatment[subject] # time coefficient g.0.true <- as.numeric( HP['g.0.true'] ) # treatment coefficient g.1.true <- -as.numeric(t)*g.0.true # intercept mu.a.true <- as.numeric( HP['mu.a.true'] ) # fixed effects b.true <- (g.0.true + g.1.true*treatment) # random effects sigma.y.true <- as.numeric( HP['sigma.y.true'] ) # residual std dev sigma.a.true <- as.numeric( HP['sigma.a.true'] ) # site std dev sigma.a0.true <- as.numeric( HP['sigma.a0.true'] ) # site:person std dev a0.true <- rnorm(J*K, 0, sigma.a0.true) a.true <- rnorm(J*K, mu.a.true + a0.true, sigma.a.true) y <- rnorm(J*K*L, a.true[subject] + b.true[subject]*time, sigma.y.true) return(data.frame( y, time, subject, treatment1, site1 )) I then fit models of the form: y ~ time + time:treatment1 + (1 | site1/subject) To the fake data. However, this approach can (but not always) produce a 'site' standard deviation approximately a factor of 10 less than in the real data. My question is - is my simulation function correct? [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.