Problems with a simulation study of the repeated experiments explanation of a 95% Confidence interval - where am I going wrong?

Question

I'm trying to write an R script to simulate the repeated experiments interpretation of a 95% confidence interval. I've found that it overestimates the proportion of times in which the true population value of a proportion is contained within the sample's 95% CI. Not a big difference - about 96% vs 95% but this interested me nonetheless.

My function takes a sample samp_n from a Bernouilli distribution with probability pop_p, and then calculates a 95% confidence interval with prop.test() using continuity correction, or more exactly with binom.test(). It returns 1 if the true population proportion pop_p is contained within the 95% CI. I've written two functions, one which uses prop.test() and one which uses binom.test() and have had similar results with both:

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

I've found that when you repeat the expriment a few thousand times, the proportion of times when the pop_p is within the 95% CI of the sample is closer to 0.96 rather than 0.95.

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

My thoughts so far about why this may be the case are

• my code is wrong (but I've checked it a lot)
• I initially thought that this was due to the normal approximation issue, but then found binom.test()

Any suggestions?

Answer 1

You're not going wrong. It simply isn't possible to construct a confidence interval for a binomial proportion which always has coverage of exactly 95% due to the discrete nature of the outcome. The Clopper-Pearson ('exact') interval is guaranteed to have coverage of at least 95%. Other intervals have coverage closer to 95% on average, when averaged over the true proportion.

I tend to favour the Jeffreys interval myself, as it has coverage close to 95% on average and (unlike the Wilson score interval) approximately equal coverage in both tails.

With only a small change in the code in the question, we can compute the exact coverage without simulation.

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

This yields the following output.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087