Mengapa interval kredibel Bayesian dalam regresi polinomial ini bias sedangkan interval kepercayaannya benar?

Pertimbangkan plot di bawah ini di mana saya mensimulasikan data sebagai berikut. Kami melihat hasil biner yang kemungkinan benar menjadi 1 ditunjukkan dengan garis hitam. Hubungan fungsional antara kovariat x dan p ( y o b s = 1 | x ) adalah polinomial urutan ke-3 dengan tautan logistik (sehingga bersifat non-linear dalam dua arah).$y_{obs}$$x$$p(y_{obs}=1 | x)$

Garis hijau adalah kecocokan regresi logistik GLM di mana diperkenalkan sebagai polinomial orde ketiga. Garis hijau putus-putus adalah interval kepercayaan 95% sekitar prediksi p ( y o b s = 1 | x , β ) , di mana β koefisien regresi dipasang. Saya menggunakan dan untuk ini.$x$$p(y_{obs}=1 | x, \hat{\beta})$$\hat{\beta}$R glmpredict.glm

Demikian pula, garis pruple adalah rata-rata posterior dengan interval kredibel 95% untuk dari model regresi logistik Bayesian menggunakan seragam sebelumnya. Saya menggunakan paket dengan fungsi untuk ini (pengaturan memberikan seragam sebelumnya tidak informatif).$p(y_{obs}=1 | x, \beta)$MCMCpackMCMClogitB0=0

Titik-titik merah menunjukkan pengamatan dalam set data yang , titik-titik hitam adalah pengamatan dengan y o b s = 0 . Perhatikan bahwa seperti biasa dalam klasifikasi / analisis diskrit y tetapi tidak p ( y o b s = 1 | x ) diamati.$y_{obs}=1$$y_{obs}=0$$y$$p(y_{obs}=1 | x)$

Beberapa hal dapat dilihat:

1. Saya disimulasikan dengan sengaja bahwa jarang di tangan kiri. Saya ingin agar interval kepercayaan dan kredibilitas melebar di sini karena kurangnya informasi (pengamatan).$x$
2. $y_{obs}=1$
3. Interval kepercayaan semakin lebar seperti yang diharapkan, sedangkan interval yang kredibel tidak . Bahkan interval kepercayaan mencakup ruang parameter yang lengkap, sebagaimana seharusnya karena kurangnya informasi.

$x$

1. Apa alasannya?
2. Langkah apa yang bisa saya ambil untuk mencapai interval kredibel yang lebih baik? (Yaitu, yang menutupi setidaknya bentuk fungsional yang sebenarnya, atau lebih baik selebar interval kepercayaan)

Kode untuk mendapatkan interval prediksi dalam grafik dicetak di sini:

fit <- glm(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
x_pred <- seq(0, 1, by=0.01)
pred <- predict(fit, newdata = data.frame(x=x_pred), se.fit = T)
plot(plogis(pred$fit), type='l')
matlines(plogis(pred$fit + pred$se.fit %o% c(-1.96,1.96)), type='l', col='black', lty=2)


library(MCMCpack)
mcmcfit <- MCMClogit(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
gibbs_samps <- as.mcmc(mcmcfit)
x_pred_dm <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=x_pred))
gibbs_preds <- apply(gibbs_samps, 1, `%*%`, t(x_pred_dm))
gibbs_pis <- plogis(apply(gibbs_preds, 1, quantile, c(0.025, 0.975)))
matlines(t(gibbs_pis), col='red', lty=2)

Akses data : https://pastebin.com/1H2iXiew terima kasih @DeltaIV dan @AdamO

$X$$X$

GLM sering binomial tidak berbeda dari GLM dengan tautan identitas kecuali bahwa varians sebanding dengan rata-rata.

$X\rightarrow -\infty$$X\rightarrow \infty$

Untuk prediksi yang sering terjadi, peningkatan proporsional deviasi kuadrat (leverage) dalam varian prediksi mendominasi kecenderungan ini. Inilah sebabnya mengapa tingkat konvergensi ke interval prediksi kira-kira sama dengan [0, 1] lebih cepat daripada konvergensi logit polinomial urutan ketiga dengan probabilitas 0 atau 1 secara tunggal.

Ini tidak demikian untuk kuantil dipasang posterior Bayesian. Tidak ada penggunaan deviasi kuadrat secara eksplisit, jadi kami hanya mengandalkan proporsi kecenderungan dominan 0 atau 1 untuk membangun interval prediksi jangka panjang.

$X$

Menggunakan kode yang saya berikan di atas kita dapatkan:

> x_pred_dom <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=c(1000)))
> gibbs_preds <- plogis(apply(gibbs_samps[1000:10000, ], 1, `%*%`, t(x_pred_dom))) # a bunch of 0/1s basically past machine precision
> prop.table(table(gibbs_preds))
gibbs_preds
         0          1 
0.97733585 0.02266415 
> 

Jadi 97,75% dari waktu, istilah polinomial ketiga adalah negatif. Ini dapat diverifikasi dari sampel Gibbs:

> prop.table(table(gibbs_samps[, 4]< 0))

 FALSE   TRUE 
0.0225 0.9775 

$X$

Di sisi lain, kesesuaian frequentist hingga 0,1 seperti yang diharapkan:

freq <- predict(fit, newdata = data.frame(x=1000), se.fit=T)
plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)

memberi:

> plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)
     [,1]
[1,]    0
[2,]    1