CODEWITHSUNDEEP
rsimulationlme4mixed-models
blazej
blazej

Reputation: 1788

Simulate data for mixed-effects model with predefined parameter

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:

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.

Upvotes: 3

Views: 2365

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226871

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

Upvotes: 3

Related Questions