How to build a random intercept and tilt in a mixed model with multiple predictors?

Question

Is it possible to construct a random interception or slope of a mixed model when it has more than one predictor?

With one predictor, I would do the following:

#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 such a model ??

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

Or using lmer

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

Should I consider all the coefficients or only those of the variables that I draw?

thanks

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 try to match a random effect for a grouping variable with only two levels - this will almost always lead to an estimate of the variance of the random zero effect, which in turn will put your predicted lines directly on top of each other - I switch 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 which makes this a little easier ...

• Predict for the range pred1 with pred2 equal to its average value, and vice versa. This is all a little smarter than it needs to be, as it generates all the values ​​for both focal predictors and speaks to ggplot at a time ...

()

 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 only on pred1 and generating predictions for the (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 a presentation, you can replace the colour aesthetics with group if all you want to do is to distinguish between groups (i.e. split separate lines), rather than identify them ...