r - 在 R 中模拟 lmer() 模型中的交互效果

Question

是否有具有以下功能的 R 包：

(1) 模拟交互变量的不同值，(2) 绘制图表，展示交互中项的不同值对 Y 的影响，以及 (3) 与拟合 lmer( ) lme4 包的功能？

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

Answer 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"))

Answer 2

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

Answer 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()

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