Masalah dengan studi simulasi penjelasan percobaan berulang tentang interval kepercayaan 95% - di mana saya salah?

Saya mencoba menulis skrip R untuk mensimulasikan interpretasi eksperimen berulang dari interval kepercayaan 95%. Saya telah menemukan bahwa itu melebih-lebihkan proporsi waktu di mana nilai populasi sebenarnya dari suatu proporsi terkandung dalam 95% CI sampel. Bukan perbedaan besar - sekitar 96% vs 95% tetapi ini tetap menarik bagi saya.

Fungsi saya mengambil sampel samp_ndari distribusi Bernoulli dengan probabilitas pop_p, dan kemudian menghitung interval kepercayaan 95% dengan prop.test()menggunakan koreksi kontinuitas, atau lebih tepatnya dengan binom.test(). Ini mengembalikan 1 jika proporsi populasi sebenarnya pop_pterkandung dalam CI 95%. Saya telah menulis dua fungsi, yang menggunakan prop.test()dan yang menggunakan binom.test()dan memiliki hasil yang serupa dengan keduanya:

in_conf_int_normal <- function(pop_p = 0.3, samp_n = 1000, correct = T){
    ## uses normal approximation to calculate confidence interval
    ## returns 1 if the CI contain the pop proportion
    ## returns 0 otherwise
    samp <- rbinom(samp_n, 1, pop_p)
    pt_result <- prop.test(length(which(samp == 1)), samp_n)
    lb <- pt_result$conf.int[1]
        ub <- pt_result$conf.int[2]
    if(pop_p < ub & pop_p > lb){
        return(1)
    } else {
    return(0)
    }
}
in_conf_int_binom <- function(pop_p = 0.3, samp_n = 1000, correct = T){
    ## uses Clopper and Pearson method
    ## returns 1 if the CI contain the pop proportion
    ## returns 0 otherwise
    samp <- rbinom(samp_n, 1, pop_p)
    pt_result <- binom.test(length(which(samp == 1)), samp_n)
    lb <- pt_result$conf.int[1]
        ub <- pt_result$conf.int[2] 
    if(pop_p < ub & pop_p > lb){
        return(1)
    } else {
    return(0)
    }
 }

Saya telah menemukan bahwa ketika Anda mengulangi percobaan beberapa ribu kali, proporsi kali ketika pop_pberada dalam CI 95% dari sampel lebih dekat ke 0,96 daripada 0,95.

set.seed(1234)
times = 10000
results <- replicate(times, in_conf_int_binom())
sum(results) / times
[1] 0.9562

Pikiranku sejauh ini tentang mengapa ini mungkin terjadi

• kode saya salah (tapi saya sudah sering memeriksanya)
• Saya awalnya berpikir bahwa ini adalah karena masalah perkiraan yang normal, tetapi kemudian ditemukan binom.test()

Ada saran?

Anda tidak salah. Sangat tidak mungkin untuk membangun interval kepercayaan untuk proporsi binomial yang selalu memiliki cakupan tepat 95% karena sifat diskrit dari hasil. Interval Clopper-Pearson ('tepat') dijamin memiliki cakupan setidaknya 95%. Interval lain memiliki cakupan mendekati rata-rata 95% , ketika dirata-ratakan atas proporsi sebenarnya.

Saya cenderung menyukai interval Jeffreys sendiri, karena memiliki cakupan mendekati rata-rata 95% dan (tidak seperti interval skor Wilson) kira-kira cakupan yang sama di kedua ekor.

Dengan hanya perubahan kecil pada kode dalam pertanyaan, kita dapat menghitung cakupan yang tepat tanpa simulasi.

p <- 0.3
n <- 1000

# Normal test
CI <- sapply(0:n, function(m) prop.test(m,n)$conf.int[1:2])
caught.you <- which(CI[1,] <= p & p <= CI[2,])
coverage.pr <- sum(dbinom(caught.you - 1, n, p))

# Clopper-Pearson
CI <- sapply(0:n, function(m) binom.test(m,n)$conf.int[1:2])
caught.you.again <- which(CI[1,] <= p & p <= CI[2,])
coverage.cp <- sum(dbinom(caught.you.again - 1, n, p))

Ini menghasilkan output berikut.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087