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
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.merMod
dokumentasi, 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 bootMer
untuk tugas ini," tetapi dengan ampun, saya tidak bisa melihat cara menggunakan bootMer
untuk melakukan ini. Tampaknya bootMer
akan 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 sec
rata-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.0
di 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 bootMer
penggunaan use.u=FALSE
secara 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 arm
paket.
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 lattice
paket:
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()
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?
predictInterval
termasuk kesalahan / ketidakpastian untuk ketentuan efek tetap dan acak. Di dalamdotplot
Anda 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 tetapfishWt
dan 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 dalampredictInterval
hasil.Jawaban:
Pertanyaan ini dan pertukaran yang sangat baik adalah dorongan untuk membuat
predictInterval
fungsi dalammerTools
paket.bootMer
adalah 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,
predictInterval
dirancang untuk menggunakanarm::sim
fungsi untuk menghasilkan distribusi parameter dalam model dan kemudian menggunakan distribusi tersebut untuk menghasilkan nilai simulasi dari respons yang diberikannewdata
oleh pengguna. Mudah digunakan - yang perlu Anda lakukan adalah:Anda dapat menentukan seluruh host nilai-nilai lain untuk
predictInterval
memasukkan 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
theta
parameter dalamlmer
objek tidak termasuk, tetapi semua variasi lain ditangkap melalui metode ini, memberikan perkiraan yang cukup baik.sumber
predictInterval()
suka efek acak bersarang? Misalnya, menggunakanmsleep
dataset dariggplot2
paket: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
devtools::install_github("jknowles/merTools")
terlebih dahulu.Lakukan ini dengan membuat
bootMer
himpunan prediksi untuk setiap replikasi bootstrap parametrik:Keluaran
bootMer
dalam objek yang tidak terlalu transparan"boot"
, tetapi kita bisa mendapatkan prediksi mentah dari$t
komponen.Seberapa sering Fish E mengalahkan Fish D?
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:
sumber
sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10
. Saat saya gunakanpredict()
, 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?use.u=TRUE
seperti pada:bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101,use.u=TRUE)
tampaknya memberi saya apa yang saya inginkan. Terima kasih!use.u
argumen untukbootMer
. 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 menggunakanuse.u=TRUE
, tetapi saya tidak berpikir itu akan selalu melakukan apa yang Anda inginkan ...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.bootMer
mengasumsikan tidak ada kesalahan dalam estimasi titik RE kami. Jika saya menggunakanuse.u=FALSE
, apakahbootMer
memperhitungkan estimasi titik RE sama sekali? TampaknyabootMer
hasil ketika menggunakanuse.u=FALSE
setara (atau, asimptotik equiv) untuk digunakanre.form=NA
dalampredict()
pernyataan. Benarkah?c(attr(ranef(lme1,condVar=TRUE)[[1]],"postVar"))
(semuanya identik dalam contoh ini), lalu sampel nilai-nilai tersebut.