在 R 中模拟 lmer() 模型中的交互效果答案

【问题标题】：Simulating an interaction effect in a lmer() model in R在 R 中模拟 lmer() 模型中的交互效果
【发布时间】：2013-04-03 06:21:26
【问题描述】：

是否有 R 包的功能可以：

(1) 模拟交互变量的不同值， (2) 绘制一个图表，说明交互作用中项的不同值对 Y 的交互作用的影响，以及 (3) 与装有 lme4 包的 lmer() 函数的模型一起工作吗？

我查看了 arm、ez、coefplot2 和 fanovaGraph 包，但找不到我想要的。

【问题讨论】：

标签： r simulation interaction lme4 lmer

【解决方案1】：

我不确定是否有包，但您可以模拟交互中不同术语的数据，然后绘制图表。这是一个通过波（即纵向）交互和绘制语法进行处理的示例。我认为这个例子背后的故事是一种提高学龄儿童口语阅读流畅度的治疗方法。通过改变 bX 的函数值来修改交互项。

library(arm)
sim1 <- function (b0=50, bGrowth=4.672,bX=15, b01=.770413, b11=.005, Vint=771, Vslope=2.24,  Verror=40.34) {

  #observation ID
  oID<-rep(1:231)
  #participant ID
  ID<-rep(1:77, each=3)
  tmp2<-sample(0:1,77,replace=TRUE,prob=c(.5,.5))
  ITT<-tmp2[ID]
  #longitudinal wave: for example 0, 4, and 7 months after treatment
  wave <-rep(c(0,4,7), 77)

  bvaset<-rnorm(77, 0, 11.58)
  bva<-bvaset[ID]

  #random effect intercept
  S.in <- rnorm(77, 0, sqrt(Vint))
  #random effect for slope
  S.sl<-rnorm(77, 0, sqrt(Vslope))
  #observation level error
  eps <- rnorm(3*77, 0, sqrt(Verror))
  #Create Outcome as product of specified model
  ORFset <- b0 + b01*bva+ bGrowth*wave +bX*ITT*wave+ S.in[ID]+S.sl[ID]*wave+eps[oID]
  #if else statement to elimiante ORF values below 0
  ORF<-ifelse(ORFset<0,0,ORFset)
  #Put into a data frame
  mydata <- data.frame( oID,ID,ITT, wave,ORF,bva,S.in[ID],S.sl[ID],eps)
    #run the model
  fit1<-lmer(ORF~1+wave+ITT+wave:ITT+(1+wave|ID),data=mydata)
  fit1

  #grab variance components
  vc<-VarCorr(fit1)

  #Select Tau and Sigma to select in the out object
  varcomps=c(unlist(lapply(vc,diag)),attr(vc,"sc")^2)
  #Produce object to output
  out<-c(coef(summary(fit1))[4,"t value"],coef(summary(fit1))[4,"Estimate"],as.numeric(varcomps[2]),varcomps[3])
  #outputs T Value, Estimate of Effect, Tau, Sigma Squared
  out
  mydata
}

mydata<-sim1(b0=50, bGrowth=4.672, bX=1.25, b01=.770413, b11=.005, Vint=771, Vslope=2.24,  Verror=40.34)
xyplot(ORF~wave,groups=interaction(ITT),data=mydata,type=c("a","p","g"))

【讨论】：

【解决方案2】：

尝试 languageR 包或效果包中的 plotLMER.fnc()。

【讨论】：

【解决方案3】：

merTools 包有一些功能可以简化此操作，但它仅适用于使用 lmer 和 glmer 对象。您可以这样做：

library(merTools)
# fit an interaction model
m1 <- lmer(y ~ studage * service + (1|d) + (1|s), data = InstEval)
# select an average observation from the model frame
examp <- draw(m1, "average")
# create a modified data.frame by changing one value
simCase <- wiggle(examp, var = "service", values = c(0, 1))
# modify again for the studage variable
simCase <- wiggle(simCase, var = "studage", values = c(2, 4, 6, 8))

在此之后，我们的模拟数据如下所示：

simCase
     y     studage service   d   s
1 3.205745       2       0 761 564
2 3.205745       2       1 761 564
3 3.205745       4       0 761 564
4 3.205745       4       1 761 564
5 3.205745       6       0 761 564
6 3.205745       6       1 761 564
7 3.205745       8       0 761 564
8 3.205745       8       1 761 564

接下来，我们需要生成预测区间，我们可以使用 merTools::predictInterval（或者没有区间，您可以使用 lme4::predict）

preds <- predictInterval(m1, level = 0.9, newdata = simCase)

现在我们得到一个 preds 对象，它是一个 3 列的 data.frame：

preds
       fit       lwr      upr
1 3.312390 1.2948130 5.251558
2 3.263301 1.1996693 5.362962
3 3.412936 1.3096006 5.244776
4 3.027135 1.1138965 4.972449
5 3.263416 0.6324732 5.257844
6 3.370330 0.9802323 5.073362
7 3.410260 1.3721760 5.280458
8 2.947482 1.3958538 5.136692

然后我们可以把它们放在一起来绘制：

library(ggplot2)
plotdf <- cbind(simCase, preds)
ggplot(plotdf, aes(x = service, y = fit, ymin = lwr, ymax = upr)) + 
  geom_pointrange() + facet_wrap(~studage) + theme_bw()

不幸的是，这里的数据导致了一个相当无趣但易于解释的情节。

【讨论】：