How to plot random intercept and slope in a mixed model with multiple predictors?

Question

Is it possible to plot the random intercept or slope of a mixed model when it has more than one predictor?

With one predictor I would do like this:

#generate one response, two predictors and one factor (random effect)
resp<-runif(100,1, 100)
pred1<-c(resp[1:50]+rnorm(50, -10, 10),resp[1:50]+rnorm(50, 20, 5))
pred2<-resp+rnorm(100, -10, 10)
RF1<-gl(2, 50)

#gamm
library(mgcv)
mod<-gamm(resp ~ pred1, random=list(RF1=~1))
plot(pred1, resp, type="n")
for (i in ranef(mod$lme)[[1]]) {
abline(fixef(mod$lme)[1]+i, fixef(mod$lme)[2])
}

#lmer
library(lme4)
mod<-lmer(resp ~ pred1 + (1|RF1))
plot(pred1, resp, type="n")
for (i in ranef(mod)[[1]][,1]) {
abline(fixef(mod)[1]+i, fixef(mod)[2])
}

But what if I have a model like this instead?:

mod<-gamm(resp ~ pred1 + pred2, random=list(RF1=~1))

Or with lmer

mod<-lmer(resp ~ pred1 + pred2 + (1|RF1))

Should I consider all the coefficients or only the ones of the variable that I'm plotting?

Thanks

Balneal

Basically, you have to decide what you want to do about the other variables. The most common procedure is to pick a reference value for one variable (e.g. pred2 equal to its mean) and plot the slope with respect to pred1 for that value. Or you could pick several values of pred2 and plot a (set of) lines for each one, possibly in separate subplots, or (ugliest) do 3D plots and plot planes resp~f(pred1,pred2) instead. – Balaklava

Thank you Ben, Sorry but I'm not sure to follow you, what do you mean exactly for "pick a reference value for one variable"? How would you do it in practice? – Balneal

Answer 1

## generate one response, two predictors and one factor (random effect)
set.seed(101)
resp <- runif(100,1,100)
pred1<- rnorm(100, 
           mean=rep(resp[1:50],2)+rep(c(-10,20),each=50),
           sd=rep(c(10,5),each=50))
pred2<- rnorm(100, resp-10, 10)

NOTE that you should probably not be trying to fit a random effect for an grouping variable with only two levels -- this will almost invariably result in an estimated random-effect variance of zero, which will in turn put your predicted lines right on top of each other -- I'm switching from gl(2,50) to gl(10,10) ...

RF1<-gl(10,10)
d <- data.frame(resp,pred1,pred2,RF1)

#lmer
library(lme4)
mod <- lmer(resp ~ pred1 + pred2 + (1|RF1),data=d)

The development version of lme4 has a predict() function that makes this a little easier ...

• Predict for a range of pred1 with pred2 equal to its mean, and vice versa. This is all a little bit cleverer than it needs to be, since it generates all the values for both focal predictors and plots them with ggplot in one go ...

()

nd <- with(d,
           rbind(data.frame(expand.grid(RF1=levels(RF1),
                      pred1=seq(min(pred1),max(pred1),length=51)),
                      pred2=mean(pred2),focus="pred1"),
                 data.frame(expand.grid(RF1=levels(RF1),
                      pred2=seq(min(pred2),max(pred2),length=51)),
                      pred1=mean(pred1),focus="pred2")))
nd$val <- with(nd,pred1[focus=="pred1"],pred2[focus=="pred2"])
pframe <- data.frame(nd,resp=predict(mod,newdata=nd))
library(ggplot2)
ggplot(pframe,aes(x=val,y=resp,colour=RF1))+geom_line()+
         facet_wrap(~focus,scale="free")

• Alternatively, focusing just on pred1 and generating predictions for a (small/discrete) range of pred2 values ...

()

nd <- with(d,
           data.frame(expand.grid(RF1=levels(RF1),
                      pred1=seq(min(pred1),max(pred1),length=51),
                      pred2=seq(-20,100,by=40))))
pframe <- data.frame(nd,resp=predict(mod,newdata=nd))
ggplot(pframe,aes(x=pred1,y=resp,colour=RF1))+geom_line()+
         facet_wrap(~pred2,nrow=1)

You might want to set scale="free" in the last facet_wrap() ... or use facet_grid(~pred2,labeller=label_both)

For presentation you might want to replace the colour aesthetic, with group, if all you want to do is distinguish among groups (i.e. plot separate lines) rather than identify them ...

Balaklava