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.

the only part of this that seems hard is specifying a Cohen's d in advance. If we ignore the random effects, can you work backward from the specified standard deviations analytically to get the required fixed effect parameters (or vice versa)? Or are you willing to approximate by this by the standard deviation-scaled effect size? — Ben Bolker
I'm not sure if I understood your question. If it's about my ability to write code or understand formulas I'm not an expert but I successfully simulated data for y ~ 1 + (1 | subject) where the ratio of random variances was as I intended (for example 50/50 or some other). I didn't know how to include different level variables and keep the random variance ratio, how to include a correlation between fixed effects and how to include Cohen's d (but I can live without the ES if you say it's hard to get) — blazej

Ben Bolker Ben Bolker · Accepted Answer · 2018-08-21T00:17:08

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

Simulate data for mixed-effects model with predefined parameter

1 Answers