Problemy z badaniem symulacyjnym powtórzonych eksperymentów z wyjaśnieniem 95% przedziału ufności - gdzie się mylę?

Próbuję napisać skrypt R, aby zasymulować interpretację powtarzanych eksperymentów z 95% przedziałem ufności. Przekonałem się, że przecenia on odsetek przypadków, w których prawdziwa wartość populacyjna części jest zawarta w 95% CI próby. Niewielka różnica - około 96% vs 95%, ale to mnie jednak zainteresowało.

Moja funkcja pobiera próbkę samp_nz rozkładu Bernoulliego z prawdopodobieństwem pop_p, a następnie oblicza 95% przedział ufności za prop.test()pomocą korekcji ciągłości, a dokładniej za pomocą binom.test(). Zwraca 1, jeśli prawdziwy odsetek populacji pop_pzawiera się w 95% CI. Napisałem dwie funkcje, jedną, która używa prop.test()i jedną, która używa binom.test()i przyniosła podobne wyniki dla obu:

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)
    }
 }

Przekonałem się, że kiedy powtórzysz eksperyment kilka tysięcy razy, odsetek przypadków, gdy pop_pmieści się w 95% CI próbki, jest bliższy 0,96 niż 0,95.

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

Moje dotychczasowe przemyślenia na temat tego, dlaczego tak się dzieje

• mój kod jest nieprawidłowy (ale często go sprawdzałem)
• Początkowo myślałem, że to z powodu normalnego problemu z przybliżeniem, ale potem znalazłem binom.test()

Jakieś sugestie?

Nie pomylisz się. Po prostu niemożliwe jest zbudowanie przedziału ufności dla proporcji dwumianowej, która zawsze ma dokładnie 95% pokrycie z powodu dyskretnej natury wyniku. Gwarantowany interwał Cloppera-Pearsona („dokładny”) wynosi co najmniej 95%. Inne przedziały mają zasięg zbliżony średnio do 95% , gdy uśrednia się je względem rzeczywistej proporcji.

Ja sam preferuję interwał Jeffreysa, ponieważ ma on zasięg średnio blisko 95% i (w przeciwieństwie do przedziału punktacji Wilsona) w przybliżeniu jednakowy zasięg w obu ogonach.

Przy niewielkiej zmianie kodu w pytaniu możemy obliczyć dokładny zasięg bez symulacji.

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))

Daje to następujące dane wyjściowe.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087