[R] Simulating data with nested random effects

[Peter Wijeratne]

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]]