Simulate data for mixed-effects model with predefined parameter

Question

I'm trying to simulate data for a model expressed with the following formula:
lme4::lmer(y ~ a + b + (1|subject), data) but with a set of given parameters:

• a <- rnorm() measured at subject level (e.g nSubjects = 50)
• y is measured at the observation level (e.g. nObs = 7 for each subject
• b <- rnorm() measured at observation level and correlated at a given r with a
• variance ratio of the random effects in lmer(y ~ 1 + (1 | subject), data) is fixed at for example 50/50 or 10/90 (and so on)
• some random noise is present (so that a full model does not explain all the variance)
• effect size of the fixed effects can be set at a predefined level (e.g. dCohen=0.5)

I played with various packages like: powerlmm, simstudy or simr but still fail to find a working solution that will accommodate the amount of parameters I'd like to define beforehand.

Also for my learning purposes I'd prefer a base R method than a package solution.

The closest example I found is a blog post by Ben Ogorek "Hierarchical linear models and lmer" which looks great but I can't figure out how to control for parameters listed above.

Any help would be appreciated. Also if there a package that I don't know of, that can do these type of simulations please let me know.

Answer 1

Some questions about the model definition:

• How do we specify a correlation between two random vectors that are different lengths? I'm not sure: I'll sample 350 values (nObs*nSubject) and throw away most of the values for the subject-level effect.
• Not sure about "variance ratio" here. By definition, the theta parameters (standard deviations of the random effects) are scaled by the residual standard deviation (sigma), e.g. if sigma=2, theta=2, then the residual std dev is 2 and the among-subject std dev is 4

Define parameter/experimental design values:

nSubjects <- 50
nObs <- 7
## means of a,b are 0 without loss of generality
sdvec <- c(a=1,b=1)
rho <- 0.5  ## correlation
betavec <- c(intercept=0,a=1,b=2)
beta_sc <- betavec[-1]*sdvec  ## scale parameter values by sd
theta <- 0.4  ## = 20/50
sigma <- 1

Set up data frame:

library(lme4)      
set.seed(101)
## generate a, b variables
mm <- MASS::mvrnorm(nSubjects*nObs,
          mu=c(0,0),
          Sigma=matrix(c(1,rho,rho,1),2,2)*outer(sdvec,sdvec))
subj <- factor(rep(seq(nSubjects),each=nObs))  ## or ?gl
## sample every nObs'th value of a
avec <- mm[seq(1,nObs*nSubjects,by=nObs),"a"]
avec <- rep(avec,each=nObs)  ## replicate
bvec <- mm[,"b"]
dd <- data.frame(a=avec,b=bvec,Subject=subj)

Simulate:

dd$y <- simulate(~a+b+(1|Subject),
               newdata=dd,
               newparams=list(beta=beta_sc,theta=theta,sigma=1),
               family=gaussian)[[1]]