[R] Simulating data with nested random effects
Peter Wijeratne
pwijer@tne @ending from gm@il@com
Wed Oct 17 21:48:25 CEST 2018
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]]
More information about the R-help
mailing list