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

[Thierry Onkelinx]

Dear Pete,

I rewrote your code to make it shorter and IMHO more readable. Meaningful
variables names require less comments on what they are. The model yields
sensible estimates on data simulated with the code below.

library(lme4)
library(dplyr)
create_fake <- function(
  n_site = 100, n_subject_site = 10, n_time = 10,
  intercept = 10, trend = 0.1, effect = 0.5,
  sigma_site = 5, sigma_subject = 2, sigma_noise = 1
){
  re.site <- rnorm(n_site, mean = 0, sd = sigma_site)
  re.subject <- rnorm(n_site * n_subject_site, mean = 0, sd = sigma_subject)

  expand.grid(
    time = seq(0, 2, length = n_time),
    site = seq_len(n_site),
    subject = seq_len(n_subject_site)
  ) %>%
    mutate(
      subject = interaction(site, subject),
      treatment = sample(0:1, size = n_site * n_subject_site,, replace =
TRUE)[subject],
      fixed = intercept + effect * treatment + trend * time,
      random = re.site[site] + re.subject[subject],
      mu = fixed + random,
      y = rnorm(n(), mean = mu, sd = sigma_noise)
    )
}
dataset <- create_fake()
m <- lmer(y ~ treatment + time + (1|site/subject), data = dataset)
summary(m)

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx using inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

<https://www.inbo.be>


Op do 18 okt. 2018 om 10:28 schreef Wijeratne, Peter <p.wijeratne using ucl.ac.uk
>:

> 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]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]