R help - Oct 2018 - Simulating data with nested random effects

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

This would almost certainly fit better on the r-sig-mixed-models
list,rather than here. You are more likely to get authoritative responses
about this specialized statistical topic there.

Also -- these are **plain text** mailing lists. Please do not post in html.

Cheers,
Bert


Bert Gunter

"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Wed, Oct 17, 2018 at 5:45 PM Peter Wijeratne <pwijeratne at gmail.com>
wrote:
> 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 at 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.
>
	[[alternative HTML version deleted]]