Как построить график случайного пересечения и наклона в смешанной модели с несколькими предикторами?

Можно ли построить случайную точку пересечения или наклон смешанной модели, когда она имеет более одного предиктора?

С одним предиктором я бы сделал так:

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

Или с лмером

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

Должен ли я учитывать все коэффициенты или только те переменные, которые я рисую?

Спасибо


r plot mixed-models random-effects lmer
person Oritteropus    schedule 14.07.2013    source источник
comment
По сути, вы должны решить, что вы хотите делать с другими переменными. Наиболее распространенная процедура состоит в том, чтобы выбрать эталонное значение для одной переменной (например, pred2, равное ее среднему значению) и построить график наклона по отношению к pred1 для этого значения. Или вы можете выбрать несколько значений pred2 и построить (набор) линий для каждого из них, возможно, в отдельных подграфиках, или (самое уродливое) вместо этого сделать трехмерные графики и построить плоскости resp~f(pred1,pred2).   -  person Ben Bolker    schedule 15.07.2013
comment
Спасибо, Бен. Извините, но я не уверен, что понимаю вас, что именно вы имеете в виду, когда выбираете ссылочное значение для одной переменной? Как бы вы поступили на практике?   -  person Oritteropus    schedule 15.07.2013

Ответы (1)


arrow_upward
5
arrow_downward 

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

Возможно, вы захотите установить scale="free" в последнем facet_wrap() ... или использовать facet_grid(~pred2,labeller=label_both)

Для презентации вы можете заменить эстетику colour на group, если все, что вы хотите сделать, это различать группы (т.е. рисовать отдельные линии), а не идентифицировать их...

person Ben Bolker    schedule 15.07.2013
comment
Очень полезно! Спасибо - person Oritteropus; 15.07.2013