CODEWITHSUNDEEP
rsimulationinteractionlme4lmer
neutral
neutral

Reputation: 127

Simulating an interaction effect in a lmer() model in R

Is there an R package with a function that can:

(1) simulate the different values of an interaction variable, (2) plot a graph that demonstrates the effect of the interaction on Y for different values of the terms in interaction, and (3) works well with the models fitted with the lmer() function of the lme4 package?

I have looked in arm, ez, coefplot2, and fanovaGraph packages, but could not find what I was looking for.

Upvotes: 1

Views: 2336

Answers (3)

jknowles
jknowles

Reputation: 479

The merTools package has some functionality to make this easier, though it only applies to working with lmer and glmer objects. Here's how you might do it:

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

After this, we have our simulated data which looks like: 

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

Next, we need to generate prediction intervals, which we can do with merTools::predictInterval (or without intervals you could use lme4::predict)

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

Now we get a preds object, which is a 3 column 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

We can then put it all together to plot:

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

Unfortunately the data here results in a rather uninteresting, but easy to interpret plot.

A faceted plot of predicted values

Upvotes: 0

user2034412
user2034412

Reputation: 4282

Try plotLMER.fnc() from the languageR package, or the effects package.

Upvotes: 1

Sam
Sam

Reputation: 381

I'm not sure about a package, but you can simulate data varying the terms in the interaction, and then graph it. Here is an example for a treatment by wave (i.e. longitudinal) interaction and the syntax to plot. I think the story behind the example is a treatment to improve oral reading fluency in school age children. The term of the interaction is modified by changing the function value for 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"))

Upvotes: 2

Related Questions