[R] Simulating data with nested random effects
Bert Gunter
bgunter@4567 @end|ng |rom gm@||@com
Thu Oct 18 04:22:03 CEST 2018
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 using 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 using 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]]
More information about the R-help
mailing list