Bagaimana cara memplot intersep dan kemiringan acak dalam model campuran dengan banyak prediktor?

Apakah mungkin untuk memplot intersep atau kemiringan acak dari model campuran jika model tersebut memiliki lebih dari satu prediktor?

Dengan satu prediktor saya akan melakukan seperti ini:

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

Namun bagaimana jika saya memiliki model seperti ini?:

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

Atau dengan lmer

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

Haruskah saya mempertimbangkan semua koefisien atau hanya variabel yang saya plot saja?

Terima kasih


r plot mixed-models random-effects lmer
person Oritteropus    schedule 14.07.2013    source sumber
comment
Pada dasarnya, Anda harus memutuskan apa yang ingin Anda lakukan terhadap variabel lainnya. Prosedur yang paling umum adalah memilih nilai referensi untuk satu variabel (misalnya pred2 sama dengan meannya) dan memplot kemiringan terhadap pred1 untuk nilai tersebut. Atau Anda dapat memilih beberapa nilai pred2 dan memplot (kumpulan) garis untuk masing-masing baris, mungkin dalam subplot terpisah, atau (yang paling jelek) membuat plot 3D dan membuat plot bidang resp~f(pred1,pred2) sebagai gantinya.   -  person Ben Bolker    schedule 15.07.2013
comment
Terima kasih Ben, Maaf tapi saya tidak yakin untuk mengikuti Anda, apa sebenarnya yang Anda maksud dengan memilih nilai referensi untuk satu variabel? Bagaimana Anda melakukannya dalam praktik?   -  person Oritteropus    schedule 15.07.2013

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

CATATAN bahwa Anda mungkin tidak mencoba menyesuaikan efek acak untuk variabel pengelompokan yang hanya memiliki dua tingkat -- hal ini hampir selalu menghasilkan perkiraan varians efek acak sebesar nol, yang pada gilirannya akan menempatkan garis prediksi Anda tepat di atas satu sama lain -- Saya beralih dari gl(2,50) ke 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)

Versi pengembangan lme4 memiliki fungsi predict() yang membuat ini sedikit lebih mudah...

  • Prediksi rentang pred1 dengan pred2 sama dengan meannya, dan sebaliknya. Ini semua sedikit lebih pintar dari yang seharusnya, karena ini menghasilkan semua nilai untuk kedua prediktor fokus dan memplotnya dengan ggplot sekaligus ...

()

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")
  • Alternatifnya, fokus hanya pada pred1 dan menghasilkan prediksi untuk rentang nilai pred2 (kecil/diskrit) ...

()

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)

Anda mungkin ingin menyetel scale="free" di facet_wrap() terakhir ... atau menggunakan facet_grid(~pred2,labeller=label_both)

Untuk presentasi Anda mungkin ingin mengganti estetika colour, dengan group, jika yang ingin Anda lakukan hanyalah membedakan antar kelompok (yaitu membuat plot baris terpisah) daripada mengidentifikasinya ...

person Ben Bolker    schedule 15.07.2013
comment
Sangat membantu! Terima kasih - person Oritteropus; 15.07.2013