Computing confidence interval of a proportion using prop.test() in R

Question

I am using prop.test() to calculate the confidence interval of a proportion. I found that prop.test() result differs from the z-score confidence interval when $\hat p$ is small, but coincide otherwise. In particular, with $\hat p=0.916$ it works well:

n = 500
p_hat = 0.916
alpha = 0.1
se = sqrt(p_hat*(1-p_hat)/n)
z = -qnorm(alpha/2)
IC = c(p_hat-z*se,p_hat+z*se); IC
## [1] 0.8955953 0.9364047

prop.test(p_hat*n,n=n,alternative="two.sided", conf.level=1-alpha, correct=FALSE)
## 0.8932886 0.9342336

But with $\hat p = 0.016$ the results differ:

IC = c(p_hat-z*se,p_hat+z*se); IC
## [1] 0.006770041 0.025229959

prop.test(p_hat*n,n=n,alternative="two.sided", conf.level=1-alpha, correct=FALSE)  
## 0.009038314 0.028171427

Do you know how in particular prop.test() calculates the confidence interval?

Answer 1

With $X = 484$ successes in $n = 500$ binomial trials, here is a comparison of the Agresti-style 95% CI with the output of 'prop.test`. The so-called Wilson CI may work even better.

n = 500;  x = 484; p0m = -1:1
x.a = x+2; n.a = n+4  # Agresti: append 2 S's & 2 F's
p.a = x.a/n.a;  se.a = sqrt(p.a*(1-p.a)/n.a)
p.a + p0m*1.96*se.a
## 0.9480839 0.9642857 0.9804876

prop.test(484, 500, correct=FALSE)$conf
## 0.9486551 0.9802085
## attr(,"conf.level")
## 0.95

Notes: (1) Getting good CIs by any method can be challenging when the proportion of successes is very near 0 or 1.

(2) The Agresti CI is not just an ad hoc trick. By conflating $1.96$ with $2$s that occur when inverting the test (and ignoring a tiny term) a 95% Agresti CI comes very close to the test-inverted CI, and yet is simple enough to remember.

(3) A frequentist re-purposing of a Bayesian probability interval based on a $\mathsf{Unif}(0,1) \equiv \mathsf{Beta}(1,1)$ prior, gives the 95% CI:

qbeta(c(.025,.975), x+1, n-x+1)
## 0.9486549 0.9801114

(4) Perhaps also see this page.