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

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 Bernoulli 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 experiment 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?

-
(+1) By the way, I re-ran your code with times=100000 a few different times and saw the same result. I'm curious to see if anyone has an explanation for this. The code is sufficiently simple that I'm pretty certain there is no coding error. Also, one run with times=1000000 gave .954931 as the result. –  Macro Jul 6 '12 at 18:00
(+1) But why do you expect to get exactly 95%? Clopper Pearson for example is guaranteed to be conservative. For your $n$ and $p$, I get that the CI should cover the true value 95.3648% of the time. –  cardinal Jul 6 '12 at 18:12
To support cardinals comment exact binomial probabilities are exact because the are based on an exact probability calculation but they do not necessarily give the exact confidence level. That is because the binomial is a discrete distribution. So Clopper-Pearson picks the endpoint for the interval so that you have the closest probability to the confidence level at or above it. This also creates a sawtoothed behavior to the power function of an exact binomial test. This odd but basic result is discussed in my paper with Christine Liu in the American Statistician (2002). –  Michael Chernick Jul 6 '12 at 18:39
Details on my paper at this link: citeulike.org/user/austin987/article/7571878 –  Michael Chernick Jul 6 '12 at 18:40
Exact Binomial CIs are "exact" because their actual performance equals their nominal performance, not because the probability calculation is "exact"! It has to be understood that a $1-\alpha$ CI must have at least a $1-\alpha$ chance of covering the true parameter, no matter what the underlying distribution is (within the assumed family). "Exact" means that the infimum of all these coverages, taken over the entire family of distributions, equals $1-\alpha$. To achieve this, the actual coverage for many of the possible distributions often must be greater than $1-\alpha$. –  whuber Jul 6 '12 at 20:57

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

-
"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" --- aside, perhaps, for the (somewhat odd) possibility of randomized intervals. (At least in that way, it can be done, though it may well be that it usually shouldn't.) –  Glen_b Aug 8 at 23:52
@Glen_b I have long been curious about the objections to randomized decisions. I believe Jack Kiefer remarked that if you're ok using randomization to collect your samples, you should have no problem using it in the decision process. If you need a decision procedure that is reproducible, documented, and difficult to cheat with, just generate any random values needed for the randomized interval before collecting the data--make it part of the design. –  whuber Aug 9 at 2:19