Interval prediksi untuk model efek gabungan lmer () di R

37

Saya ingin mendapatkan interval prediksi sekitar prediksi dari model lmer (). Saya telah menemukan beberapa diskusi tentang ini:

http://rstudio-pubs-static.s3.amazonaws.com/24365_2803ab8299934e888a60e7b16113f619.html

http://glmm.wikidot.com/faq

tetapi mereka tampaknya tidak memperhitungkan ketidakpastian efek acak.

Inilah contoh spesifik. Saya balap ikan emas. Saya memiliki data pada 100 balapan terakhir. Saya ingin memprediksi 101, dengan mempertimbangkan ketidakpastian estimasi ET saya, dan perkiraan FE. Saya termasuk intersepsi acak untuk ikan (ada 10 ikan berbeda), dan efek tetap untuk berat (ikan kurang berat lebih cepat).

library("lme4")

fish <- as.factor(rep(letters[1:10], each=100))
race <- as.factor(rep(900:999, 10))
oz <- round(1 + rnorm(1000)/10, 3)
sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10

fishDat <- data.frame(fishID = fish, 
      raceID = race, fishWt = oz, time = sec)
head(fishDat)
plot(fishDat$fishID, fishDat$time)

lme1 <- lmer(time ~ fishWt + (1 | fishID), data=fishDat)
summary(lme1)

Sekarang, untuk memprediksi balapan ke-101. Ikan-ikan telah ditimbang dan siap untuk pergi:

newDat <- data.frame(fishID = letters[1:10], 
    raceID = rep(1000, 10),
    fishWt = 1 + round(rnorm(10)/10, 3))
newDat$pred <- predict(lme1, newDat)
newDat

   fishID raceID fishWt     pred
1       a   1000  1.073 10.15348
2       b   1000  1.001 10.20107
3       c   1000  0.945 10.25978
4       d   1000  1.110 10.51753
5       e   1000  0.910 10.41511
6       f   1000  0.848 10.44547
7       g   1000  0.991 10.68678
8       h   1000  0.737 10.56929
9       i   1000  0.993 10.89564
10      j   1000  0.649 10.65480

Ikan D benar-benar membiarkan dirinya sendiri (1,11 oz) dan sebenarnya diprediksi kalah dari Ikan E dan Ikan F, keduanya lebih baik daripada sebelumnya. Namun, sekarang saya ingin mengatakan, "Ikan E (berat 0,91 oz) akan mengalahkan Ikan D (berat 1,11 oz) dengan probabilitas p." Apakah ada cara untuk membuat pernyataan seperti itu menggunakan lme4? Saya ingin probabilitas saya untuk memperhitungkan ketidakpastian saya dalam efek tetap, dan efek acak.

Terima kasih!

PS melihat predict.merModdokumentasi, itu menunjukkan "Tidak ada pilihan untuk menghitung kesalahan standar prediksi karena sulit untuk mendefinisikan metode yang efisien yang menggabungkan ketidakpastian dalam parameter varians; kami merekomendasikan bootMeruntuk tugas ini," tetapi dengan ampun, saya tidak bisa melihat cara menggunakan bootMeruntuk melakukan ini. Tampaknya bootMerakan digunakan untuk mendapatkan interval kepercayaan bootstrap untuk estimasi parameter, tapi saya bisa saja salah.

DIPERBARUI T:

OK, saya pikir saya mengajukan pertanyaan yang salah. Saya ingin dapat mengatakan, "Ikan A, dengan berat badan, akan memiliki waktu lomba yaitu (lcl, ucl) 90% dari waktu."

Dalam contoh yang saya paparkan, Ikan A, dengan berat 1,0 ons, akan memiliki waktu lomba 9 + 0.1 + 1 = 10.1 secrata-rata, dengan standar deviasi 0,1. Dengan demikian, waktu perlombaan yang diamati akan berada di antara

x <- rnorm(mean = 10.1, sd = 0.1, n=10000)
quantile(x, c(0.05,0.50,0.95))
       5%       50%       95% 
 9.938541 10.100032 10.261243 

90% dari waktu. Saya ingin fungsi prediksi yang mencoba memberi saya jawaban itu. Pengaturan semua fishWt = 1.0di newDat, kembali menjalankan sim, dan menggunakan (seperti yang disarankan oleh Ben Bolker di bawah ini)

predFun <- function(fit) {
  predict(fit,newDat)
}
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = FALSE)
predMat <- bb$t

memberi

> quantile(predMat[,1], c(0.05,0.50,0.95))
      5%      50%      95% 
10.01362 10.55646 11.05462 

Ini sepertinya berpusat pada rata-rata populasi? Seolah tidak memperhitungkan efek FishID? Saya pikir mungkin itu masalah ukuran sampel, tetapi ketika saya menabrak jumlah balapan yang diamati dari 100 menjadi 10.000, saya masih mendapatkan hasil yang serupa.

Saya akan mencatat bootMerpenggunaan use.u=FALSEsecara default. Di sisi lain, menggunakan

bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = TRUE)

memberi

> quantile(predMat[,1], c(0.05,0.50,0.95))
      5%      50%      95% 
10.09970 10.10128 10.10270 

Interval itu terlalu sempit, dan tampaknya menjadi interval kepercayaan untuk waktu rata-rata Fish A. Saya ingin interval kepercayaan untuk waktu lomba yang diamati Fish A, bukan waktu perlombaan rata-rata. Bagaimana saya bisa mendapatkannya?

UPDATE 2, HAMPIR:

Saya pikir saya menemukan apa yang saya cari di Gelman and Hill (2007) , halaman 273. Perlu memanfaatkan armpaket.

library("arm")

Untuk Ikan A:

x.tilde <- 1    #observed fishWt for new race
sigma.y.hat <- sigma.hat(lme1)$sigma$data        #get uncertainty estimate of our model
coef.hat <- as.matrix(coef(lme1)$fishID)[1,]    #get intercept (random) and fishWt (fixed) parameter estimates
y.tilde <- rnorm(1000, coef.hat %*% c(1, x.tilde), sigma.y.hat) #simulate
quantile (y.tilde, c(.05, .5, .95))

  5%       50%       95% 
 9.930695 10.100209 10.263551 

Untuk semua ikan:

x.tilde <- rep(1,10)  #assume all fish weight 1 oz
#x.tilde <- 1 + rnorm(10)/10  #alternatively, draw random weights as in original example
sigma.y.hat <- sigma.hat(lme1)$sigma$data
coef.hat <- as.matrix(coef(lme1)$fishID)
y.tilde <- matrix(rnorm(1000, coef.hat %*% matrix(c(rep(1,10), x.tilde), nrow = 2 , byrow = TRUE), sigma.y.hat), ncol = 10, byrow = TRUE)
quantile (y.tilde[,1], c(.05, .5, .95))
       5%       50%       95% 
 9.937138 10.102627 10.234616 

Sebenarnya, ini mungkin bukan yang saya inginkan. Saya hanya memperhitungkan ketidakpastian keseluruhan model. Dalam situasi di mana saya memiliki, katakanlah, 5 ras yang diamati untuk Fish K dan 1000 ras yang diamati untuk Fish L, saya pikir ketidakpastian yang terkait dengan prediksi saya untuk Fish K harus jauh lebih besar daripada ketidakpastian yang terkait dengan prediksi saya untuk Fish L.

Akan melihat lebih jauh ke Gelman dan Hill 2007. Saya merasa saya akhirnya harus beralih ke BUGS (atau Stan).

UPDATE THE 3:

Mungkin saya membuat konsep hal-hal buruk. Menggunakan predictInterval()fungsi yang diberikan oleh Jared Knowles dalam jawaban di bawah ini memberikan interval yang tidak seperti yang saya harapkan ...

library("lattice")
library("lme4")
library("ggplot2")

fish <- c(rep(letters[1:10], each = 100), rep("k", 995), rep("l", 5))
oz <- round(1 + rnorm(2000)/10, 3)
sec <- 9 + c(rep(1:10, each = 100)/10,rep(1.1, 995), rep(1.2, 5)) + oz + rnorm(2000)

fishDat <- data.frame(fishID = fish, fishWt = oz, time = sec)
dim(fishDat)
head(fishDat)
plot(fishDat$fishID, fishDat$time)

lme1 <- lmer(time ~ fishWt + (1 | fishID), data=fishDat)
summary(lme1)
dotplot(ranef(lme1, condVar = TRUE))

Saya telah menambahkan dua ikan baru. Ikan K, untuk siapa kami telah mengamati 995 ras, dan Ikan L, untuk siapa kami telah mengamati 5 ras. Kami telah mengamati 100 balapan untuk Fish AJ. Saya cocok sama lmer()seperti sebelumnya. Melihat dotplot()dari latticepaket:

Perkiraan FishID

Secara default, dotplot()mengatur ulang efek acak dengan estimasi titik mereka. Perkiraan untuk Ikan L ada di baris teratas, dan memiliki interval kepercayaan yang sangat luas. Ikan K ada di baris ketiga, dan memiliki interval kepercayaan yang sangat sempit. Ini masuk akal bagi saya. Kami memiliki banyak data tentang Fish K, tetapi tidak banyak data tentang Fish L, jadi kami lebih percaya diri dalam perkiraan kami tentang kecepatan renang Fish K yang sebenarnya. Sekarang, saya akan berpikir ini akan menyebabkan interval prediksi sempit untuk Fish K, dan interval prediksi lebar untuk Fish L saat menggunakan predictInterval(). Howeva:

newDat <- data.frame(fishID = letters[1:12],
                     fishWt = 1)

preds <- predictInterval(lme1, newdata = newDat, n.sims = 999)
preds
ggplot(aes(x=letters[1:12], y=fit, ymin=lwr, ymax=upr), data=preds) +
  geom_point() + 
  geom_linerange() +
  labs(x="Index", y="Prediction w/ 95% PI") + theme_bw()

Interval Prediksi untuk Ikan

Semua interval prediksi tersebut tampaknya memiliki lebar yang identik. Mengapa prediksi kami untuk Fish K tidak menyempit yang lain? Mengapa prediksi kami untuk Fish L tidak lebih dari yang lain?

mungkin
sumber
1
predictIntervaltermasuk kesalahan / ketidakpastian untuk ketentuan efek tetap dan acak. Di dalam dotplotAnda hanya melihat ketidakpastian karena porsi acak dari prediksi, pada dasarnya ketidakpastian di sekitar perkiraan intersep spesifik ikan. Jika model Anda memiliki banyak ketidakpastian dalam parameter tetap fishWtdan parameter ini mendorong sebagian besar nilai yang diprediksi, maka ketidakpastian di sekitar setiap intersep ikan tertentu sepele dan Anda tidak akan melihat perbedaan besar dalam lebar interval. Kita harus membuat ini lebih jelas dalam predictIntervalhasil.
Mengetahui

Jawaban:

19

Pertanyaan ini dan pertukaran yang sangat baik adalah dorongan untuk membuat predictIntervalfungsi dalam merToolspaket. bootMeradalah jalan yang harus ditempuh, tetapi untuk beberapa masalah tidak layak secara komputasional untuk menghasilkan reparasi bootstrap seluruh model (dalam kasus-kasus di mana model tersebut besar).

Dalam kasus tersebut, predictIntervaldirancang untuk menggunakan arm::simfungsi untuk menghasilkan distribusi parameter dalam model dan kemudian menggunakan distribusi tersebut untuk menghasilkan nilai simulasi dari respons yang diberikan newdataoleh pengguna. Mudah digunakan - yang perlu Anda lakukan adalah:

library(merTools)
preds <- predictInterval(lme1, newdata = newDat, n.sims = 999)

Anda dapat menentukan seluruh host nilai-nilai lain untuk predictIntervalmemasukkan pengaturan interval untuk interval prediksi, memilih apakah akan melaporkan rata-rata atau median dari distribusi, dan memilih apakah akan memasukkan varian residu dari model atau tidak.

Ini bukan interval prediksi penuh karena variabilitas thetaparameter dalam lmerobjek tidak termasuk, tetapi semua variasi lain ditangkap melalui metode ini, memberikan perkiraan yang cukup baik.

yahudi
sumber
3
Ini terlihat luar biasa! Membaca sketsa sekarang. Terima kasih!
mungkin
Interval prediksi tidak seperti yang saya harapkan. Lihat pembaruan 3 di atas.
senang hati
Tidak predictInterval()suka efek acak bersarang? Misalnya, menggunakan msleepdataset dari ggplot2paket: mod <- lmer(sleep_total ~ bodywt + (1|vore/order), data=msleep); predInt <- predictInterval(merMod=mod, newdata=msleep) Mengembalikan kesalahan:Error in '[.data.frame'(newdata, , j) : undefined columns selected
hossibley
Saya yakin itu tidak suka efek bersarang. Saya tidak berpikir kami memiliki tes apa pun di ruang tes kami. Saya akan mengajukan masalah pada GitHub untuk melihatnya. Saya juga merekomendasikan untuk mencoba versi dev dari GitHub devtools::install_github("jknowles/merTools")terlebih dahulu.
Mengetahui
2
Sebagai pembaruan, versi pengembangan terbaru dari merTools memungkinkan efek bersarang. Ini akan segera didorong ke CRAN.
jknowles
15

Lakukan ini dengan membuat bootMerhimpunan prediksi untuk setiap replikasi bootstrap parametrik:

predFun <- function(fit) {
    predict(fit,newDat)
}
bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101)

Keluaran bootMerdalam objek yang tidak terlalu transparan "boot", tetapi kita bisa mendapatkan prediksi mentah dari $tkomponen.

Seberapa sering Fish E mengalahkan Fish D?

predMat <- bb$t
dim(predMat) ## 200 rows (PB reps) x 10 (predictions)

Waktu ikan E ada di kolom 5, waktu ikan D ada di kolom 4, jadi kita hanya perlu tahu proporsi bahwa kolom 5 kurang dari kolom 4:

mean(predMat[,5]<predMat[,4])  ## 0.57
Ben Bolker
sumber
Saya mendapatkan beberapa hasil yang tidak terduga. Jika saya menetapkan fishWt = 1 untuk semua ikan di NewDat, saya akan mengharapkan rata-rata / median waktu untuk Ikan A menjadi ~ 10.1, Ikan B ~ 10.2, ..., Ikan J ~ 11.0 (karena waktu mereka dalam data pelatihan adalah didefinisikan sebagai:) sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10. Saat saya gunakan predict(), waktu prediksi untuk Ikan A, E, dan J adalah 10,09, 10,49, dan 10,99, seperti yang diharapkan. Namun, waktu rata-rata untuk metode bootMer yang Anda descirbe adalah: 10,52, 10,59, dan 10,50. Saya akan mengharapkan lebih banyak kesepakatan?
mungkin
Menggunakan use.u=TRUEseperti pada: bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101,use.u=TRUE)tampaknya memberi saya apa yang saya inginkan. Terima kasih!
mungkin
OKE, ini agak rumit. Anda perlu melihat use.uargumen untuk bootMer. Pertanyaannya adalah, ketika Anda mengatakan "ketidakpastian dalam efek tetap dan efek acak", apa yang Anda maksud dengan 'efek acak'? Apakah yang Anda maksud ketidakpastian dalam varians efek-acak , atau dalam mode kondisional (yaitu efek spesifik ikan)? Anda dapat menggunakan use.u=TRUE, tetapi saya tidak berpikir itu akan selalu melakukan apa yang Anda inginkan ...
Ben Bolker
Jika saya menggunakan use.u=TRUE, maka "nilai u [tetap] tetap pada nilai estimasi mereka". Saya menafsirkan ini sebagai makna, apa pun perkiraan titik efek acak kami untuk Ikan A, diambil sebagai Kebenaran Jujur Tuhan, jika Anda mau. bootMermengasumsikan tidak ada kesalahan dalam estimasi titik RE kami. Jika saya menggunakan use.u=FALSE, apakah bootMermemperhitungkan estimasi titik RE sama sekali? Tampaknya bootMerhasil ketika menggunakan use.u=FALSEsetara (atau, asimptotik equiv) untuk digunakan re.form=NAdalam predict()pernyataan. Benarkah?
mungkin
1
Saya pikir ini bukan ATM yang diimplementasikan, tetapi Anda dapat mengekstraksi varian bersyarat dari mode kondisional / BLUP melalui c(attr(ranef(lme1,condVar=TRUE)[[1]],"postVar"))(semuanya identik dalam contoh ini), lalu sampel nilai-nilai tersebut.
Ben Bolker