Проблемы с имитационным исследованием объяснения повторных экспериментов с 95% доверительным интервалом - где я ошибаюсь?

Я пытаюсь написать сценарий R для имитации повторяющихся экспериментов интерпретации 95% доверительного интервала. Я обнаружил, что он переоценивает долю случаев, когда истинная популяционная ценность доли содержится в 95% ДИ выборки. Не большая разница - около 96% против 95%, но это, тем не менее, меня заинтересовало.

Моя функция берет выборку samp_nиз распределения Бернулли с вероятностью pop_p, а затем вычисляет 95% доверительный интервал с prop.test()использованием коррекции непрерывности или, точнее, с binom.test(). Возвращает 1, если истинная доля населения pop_pсодержится в 95% ДИ. Я написал две функции, одна из которых использует, prop.test()а другая использует binom.test()и имеет схожие результаты с обеими:

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

Я обнаружил, что когда вы повторяете эксперимент несколько тысяч раз, доля случаев, когда pop_pДИ находится в пределах 95% ДИ образца, ближе к 0,96, чем к 0,95.

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

Мои мысли до сих пор о том, почему это может быть так

• мой код неверен (но я много проверял)
• Сначала я думал, что это связано с проблемой нормального приближения, но потом нашел binom.test()

Какие-либо предложения?

Вы не ошибетесь. Просто невозможно построить доверительный интервал для биномиальной пропорции, которая всегда имеет охват ровно 95% из-за дискретного характера результата. Интервал Клоппера-Пирсона («точный») гарантированно охватит не менее 95%. Другие интервалы охватывают ближе к 95% в среднем , когда усреднены по истинной пропорции.

Я склоняюсь к тому, чтобы отдать предпочтение интервалу Джеффриса, поскольку его охват в среднем близок к 95% и (в отличие от интервала оценок Уилсона) примерно одинаковый охват в обоих хвостах.

Только с небольшим изменением кода в вопросе, мы можем вычислить точное покрытие без моделирования.

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

Это дает следующий вывод.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087