使用预定义参数模拟混合效应模型的数据

Simulate data for mixed-effects model with predefined parameter

我正在尝试为用以下公式表示的模型模拟数据:
lme4::lmer(y ~ a + b + (1|subject), data) 但具有一组给定参数:

我尝试过各种包,例如:powerlmmsimstudysimr,但仍然无法找到适合我想要定义的参数数量的有效解决方案事先.

另外出于我的学习目的,我更喜欢基本的 R 方法而不是包解决方案。

我找到的最接近的示例是 Ben Ogorek "Hierarchical linear models and lmer" 的博客 post,它看起来不错,但我不知道如何控制上面列出的参数。

如有任何帮助,我们将不胜感激。 另外,如果有我不知道的包可以进行这些类型的模拟,请告诉我。

关于模型定义的一些问题:

  • 我们如何指定两个不同长度的随机向量之间的相关性?我不确定:我将采样 350 个值 (nObs*nSubject) 并丢弃大部分值以获得主题级效果。
  • 不确定 "variance ratio" 这里。根据定义,theta 参数(随机效应的标准差)按残差标准差 (sigma) 缩放,例如如果 sigma=2theta=2,则残差标准差为 2,主体间标准差为 4

定义parameter/experimental设计值:

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

设置数据框:

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)

模拟:

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