[R-sig-ME] Simulating data with nested random effects

[Wijeratne, Peter]

Dear r-sig-mixed-models,


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?


Note - I can generate data and provide the full code if required.


Thanks in advance for any help.

	[[alternative HTML version deleted]]