在R中模拟lmer（）模型中的交互效果

3 回答

• 0

  我不确定软件包，但您可以模拟不同交互中的条款的数据，然后绘制图形 . 以下是通过波（即纵向）相互作用和绘制语法进行治疗的示例 . 我认为这个例子背后的故事是一种提高学龄儿童口语阅读流利度的方法 . 通过更改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"))
  

  回复于 2024-05-03T19:43:29+08:00
• 2

  尝试使用languageR包中的plotLMER.fnc（）或效果包 .

  回复于 2024-05-03T19:43:29+08:00
• 1

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

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

  

  回复于 2024-05-03T19:43:29+08:00