If all tests are negative, the positive rate is zero. And its confidence interval?
Some time ago a colleague presented me a simple statistics question which, as usual, turned out to be quite interesting and intriguing. They had run 23 tests that were all negative, and they wanted to have a confidence interval for the proportion of positive outcomes. So, to put it in statistical terms, we have observations that can be 0 or 1. Each observation has a probability of being positive. In formulas
We know that the number of successes is binomially distributed, thus a sufficient statistics is and that the maximum likelihood estimator for the binomial proportion is
If the observations are all negative, the estimate for the rate is clearly zero, but what about its confidence interval? I had never been asked to estimate the confidence interval for the binomial distribution, so I was totally unprepared (shame on me!). Quite instinctively, I computed the probability to observe all 23 results negative for a given , set this to 5% and solved for . Result, 12.2%. In formula
which gives, solving for ,
In other words, what I did was to compute the highest for which I would be very surprised (5%) to see all negative outcomes. This is equivalent to being quite sure (95% sure) to observe at least one positive outcome.
My colleague had found the formula for the Clopper-Pearson interval and applied it to their case: different result. Of course, I thought: I computed the probability to observe data as extreme or more extreme than those observed for a given value of the parameter. I did a hypothesis test, not a confidence interval.
It would have been the end of the story, had I not tried this online calculator. It reports to be the 90% interval for . This was no coincidence, as the plot below shows.
In other words, for any number of tests (at least, between five and thirty) my estimate (violet points) matches the upper limit of the 90% confidence interval computed with the Clopper-Pearson method (magenta line).
library(ggplot2) library(ggthemr) ggthemr('solarized', type='outer') # computing 90% confidence interval with "exact" meaning Clopper-Pearson library(binom) cis <- binom.confint(0, 5:30, conf.level = 0.9, methods="exact") # max theta for which n all negative outcomes have probability 5% p_all_neg <- 1-0.05**(1./5:30) df <- data.frame(x=5:30, y1=cis$upper, y2=p_all_neg) p <- ggplot(df, aes(x)) + geom_line(y = df$y1, show.legend = TRUE, colour='#d33682') + geom_point(y = df$y2, show.legend = TRUE, colour='#6c71c4') + xlab("trials") + ylab(expression(theta)) + ylim(0.0, 0.5) + theme(axis.title.x = element_text(face="bold", size=18), axis.text.x = element_text(angle=0, vjust=0.0, size=14), axis.title.y = element_text(face="bold", size=18), axis.text.y = element_text(angle=0, vjust=0.0, size=14)) print(p)
Finding confidence intervals
There is a one-to-one correspondence between confidence intervals and hypothesis tests. As a matter of fact, confidence sets are found by inverting a test statistics.
Let’s start revisiting a few concepts, following the classics. We have a hypothesis for a parameter of interest $\theta$. The hypothesis says that it has a certain value
We have data and a test statistics telling us whether to reject the hypothesis or not. The acceptance region is the region of the sample space for which we do not reject at a level . In symbols
Now, for each realisation of the data , we take the values of for which the hypothesis is accepted. Then we have built a confidence set for . We define as the set in parameter
Then, is a confidence set. This follows quite immediately from the definition of acceptance region above
There are several estimators for the confidence interval of the binomial proportion. The advantage of this one is that it is exact, rather than based on the normal approximation (see the Wikipedia page linked above). The disadvantage is that it is conservative, i.e. there might be a smaller interval with the same confidence level. The estimated interval is defined as
From this definition we reconcile the fact that my estimate at 95% coincides with the Clopper-Pearson at 90%. In fact, since we are in the special case of , we can write the Clopper-Pearson as
By taking in the Clopper-Pearson estimate, we have the formula I used for my estimate.
A final observation
What is most interesting to my colleague? The hypothesis test says that those data would already exclude or higher at 5%, but the 95% interval according to Clopper-Pearson is . As we saw they are two different but intimately related things. What is more important to you depends largely on your taste. We also had the chance to underline something that is often neglected: finding an confidence interval doesn’t mean we found the smallest interval that contains the true value with probability . If you want to investigate an extreme consequence of this, you can visit the ultimate confidence intervals calculator.
Header image from Flickr