Interval prediksi untuk proporsi keberhasilan di masa depan dalam pengaturan Binomial

Misalkan saya cocok dengan regresi Binomial dan mendapatkan estimasi titik dan matriks varians-kovarians dari koefisien regresi. Itu akan memungkinkan saya untuk mendapatkan CI untuk proporsi keberhasilan yang diharapkan dalam percobaan di masa depan, , tetapi saya membutuhkan CI untuk proporsi yang diamati. Ada beberapa jawaban terkait yang diposting, termasuk simulasi (misalkan saya tidak ingin melakukan itu) dan tautan ke Krishnamoorthya et al (yang tidak cukup menjawab pertanyaan saya).$p$

Alasan saya adalah sebagai berikut: jika kita hanya menggunakan model Binomial, kita dipaksa untuk mengasumsikan bahwa diambil sampel dari distribusi Normal (dengan Wald CI yang sesuai) dan oleh karena itu tidak mungkin untuk mendapatkan CI untuk proporsi yang diamati dalam bentuk tertutup. Jika kita mengasumsikan bahwa diambil sampel dari distribusi beta, maka segala sesuatunya jauh lebih mudah karena jumlah keberhasilan akan mengikuti distribusi Beta-Binomial. Kita harus mengasumsikan bahwa tidak ada ketidakpastian dalam estimasi parameter beta, dan .$p$$p$$\alpha$$\beta$

Ada tiga pertanyaan:

1) Yang teoretis: apakah boleh menggunakan estimasi titik parameter beta saja? Saya tahu bahwa untuk membangun CI untuk pengamatan masa depan dalam regresi linier berganda

$Y = x'\beta + \epsilon, \epsilon \sim N(0, \sigma^2)$

mereka melakukan itu varian istilah kesalahan kesalahan, . Saya mengambilnya (mengoreksi saya jika saya salah) bahwa pembenarannya adalah bahwa dalam praktiknya diperkirakan dengan presisi yang jauh lebih besar daripada koefisien regresi dan kami tidak akan mendapatkan banyak dengan mencoba memasukkan ketidakpastian . Apakah justifikasi serupa berlaku untuk parameter beta yang diperkirakan, dan ?$\sigma^2$$\sigma^2$$\sigma^2$$\alpha$$\beta$

2) Paket apa yang lebih baik (R: gamlss-bb, betareg, aod?; Saya juga memiliki akses ke SAS).

3) Dengan mengestimasi parameter beta, apakah ada jalan pintas (perkiraan) untuk mendapatkan kuantil (2,5%, 97,5%) untuk menghitung keberhasilan di masa depan atau, lebih baik lagi, untuk proporsi keberhasilan di masa mendatang di bawah distribusi Beta-Binomial.

Saya akan membahas ketiga bagian pertanyaan ini.

Ada dua masalah yang tergabung, pertama adalah metode yang Anda gunakan agar sesuai dengan model regresi dalam kasus ini. Yang kedua adalah cara memperkirakan taksiran dari taksiran Anda untuk memperkirakan taksiran baru.

jika variabel respons Anda terdistribusi secara binerial, Anda biasanya akan menggunakan regresi logistik atau regresi probit (glm dengan cdf normal sebagai fungsi tautan).

Jika Anda melakukan regresi logistik, ambillah respons sebagai rasio jumlah yang diamati dibagi dengan batas atas yang diketahui, yaitu . Kemudian bawa prediktor / kovariat Anda dan masukkan ini ke dalam panggilan R ke fungsi glm. Objek yang dikembalikan memiliki semua yang Anda butuhkan untuk melakukan sisa perhitungan Anda. $y_i/n_i$

x<- rnorm(100, sd=2)
prob_true <- 1/(1+exp(-(1+5*x)))
counts <- rbinom(100, 50,prob_true)
print(d.AD <- data.frame(counts,x))
glm.D93 <- glm(counts/50 ~ x, family = binomial() )

Untuk model regresi linier , rumus untuk interval prediksi adalah:

$\hat{y}_i \pm t_{n-p}s_y\sqrt{1+\frac{1}{n}+\frac{(x_i-\bar{x})^2}{(n-1)s^2_x}}$

Anda dapat menggunakan model regresi linier sebagai perkiraan untuk glm. Untuk melakukan ini, Anda akan menggunakan rumus regresi linier untuk kombinasi linier prediktor sebelum Anda melakukan transformasi tautan terbalik untuk mendapatkan probabilitas kembali pada skala 0-1. Kode untuk melakukan ini dimasukkan ke fungsi predict.glm () R. Berikut adalah beberapa contoh kode yang juga akan membuat plot yang bagus. ( EDIT : Kode ini untuk interval kepercayaan, bukan untuk interval prediksi)

y_hat <- predict(glm.D93, type="link", se.fit=TRUE)
t_np<- qt(.975, 100-2, ncp=0)

ub <- y_hat$fit + t_np * y_hat$se.fit
lb <- y_hat$fit - t_np * y_hat$se.fit

point <- y_hat$fit

p_hat <- glm.D93$family$linkinv(point)
p_hat_lb <- glm.D93$family$linkinv(lb)
p_hat_ub <- glm.D93$family$linkinv(ub)

plot(x,p_hat)
points(x, p_hat_ub, col='red')
points(x, p_hat_lb, col='blue')

Anda dapat melakukan hal yang sama untuk glm apa pun, misalnya Poisson, invers Gaussian, gamma, dll. Dalam setiap kasus lakukan interval prediksi pada skala kombinasi linear dari prediktor. Setelah Anda mendapatkan dua titik akhir dari interval prediksi, Anda mengonversi titik akhir ini melalui tautan terbalik. Untuk masing-masing glm yang saya sebutkan tautan terbalik mungkin berbeda dari case logit yang saya tulis di sini. Semoga ini membantu.