如何在具有多个预测变量的混合模型中绘制随机截距和斜率？

Question

当混合模型拥有多个预测变量时，是否可以绘制随机截距或斜率?

如果只有一个预测变量，我会这样做:

#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])
}

但是如果我有这样的一个模型呢？：

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

使用lmer模块进行处理。

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

我应该考虑所有系数还是只考虑我正在绘制的变量的系数？

谢谢

- Oritteropus

---

基本上，你必须决定如何处理其他变量。最常见的方法是选择一个变量的参考值（例如pred2等于其平均值），并绘制相对于该值的斜率与pred1的关系图。或者您可以选择几个pred2的值，并为每个值绘制一条（组）线，可能在单独的子图中，或者（最丑陋的）进行3D绘图，并绘制平面resp〜f（pred1，pred2）。 - Ben Bolker

谢谢Ben，抱歉我不确定你的意思是什么，“为一个变量选择一个参考值”是什么意思？在实践中你会怎么做？ - Oritteropus

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)

注意，如果您尝试为仅具有两个级别的分组变量拟合随机效应，则很可能会导致估计的随机效应方差为零，从而将预测线放置在彼此重叠的位置上。我将从gl(2,50)切换到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)

lme4的开发版本有一个predict()函数，使这变得更容易...

• 针对一系列pred1进行预测，其中pred2等于其均值，反之亦然。 这比必要的聪明一些，因为它会一次性生成两个焦点预测器的所有值，并使用ggplot将它们绘制出来...

()

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

• 或者，只专注于pred1并生成对一组（小/离散）pred2值的预测...

()

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)

你可能想在最后一个facet_wrap()中设置scale="free"，或者使用facet_grid(~pred2,labeller=label_both)

如果你只是想区分组别（即绘制单独的线条），而不是识别它们，那么为了演示，你可能想用group代替colour配色美学...

- Ben Bolker