Binomial Confidence Intervals - Bayes Jeffrey's Prior vs Agresti-Coull method

The R binom library has several confidence intervals to choose from for Binomial distributions.

The Bayes method uses the Beta distribution. According to the binom documentation:

The default prior is Jeffrey's prior which is a Beta(0.5, 0.5) distribution. Thus the posterior mean is (x + 0.5)/(n + 1).

p|x ~ Beta(x + prior.shape1, n - x + prior.shape2)

The prior.shape1 and prior.shape2 can be passed in like so:

binom.bayes(x, n,
conf.level = 0.95,
type = c("highest", "central"),
prior.shape1 = 0.5,
prior.shape2 = 0.5,
tol = .Machine$double.eps^0.5, maxit = 1000, ...) Remembering that the default Bayes formula is (x + 0.5)/(n + 1), what would the shape parameters be to replicate Agresti Coull? (So instead of 0.5 and 1, they'd be replaced with 1/2*z^2 and z^2 respectively.) Question #1, what would the shape1 and shape2 parameters of the Beta distribution be to match Agresti Coull? Academically, Agresti-Coull confidence interval is considered a Bayesian method. The Agresti-Coull Interval specifies prior knowledge of z^2 for typically 3.8416 or essentially 4 given the rule of thumb "add 2 successes and 2 failures". So for this analysis, I put z=2. Notice that the Agresti-Coull more closely matches the "exact" method than does Bayes set to Jeffrey's Prior (shape1,shape2)=(1/2,1/2). If one notices the w= values, those are the necessary weight to replicate the "exact" method as per Agresti-Coull method. (w=z^2) UPDATE ON QUESTION #1 Looking at the source code for binom.bayes, one sees: a <- x + prior.shape1 b <- n - x + prior.shape2 p <- a/(a + b) Where p is the suggested mean given by the method. Expanding p, we get p <- (x + prior.shape1) / (n + prior.shape1 + prior.shape2) Notice that the x cancels out in the denominator. That implies that prior.shape1 <- prior.shape2 <- (z^2)/2 as a possibility. Given ptilde=function(x,n,z=2){ n=n+z*z p=1/n*(x+0.5*z*z) p } we get: > binom.bayes(0,25,prior.shape1 = (qnorm(0.95+.05/2)^2)/2,prior.shape2 = (qnorm(0.95+.05/2)^2)/2) method x n shape1 shape2 mean lower upper sig 1 bayes 0 25 1.920729 26.92073 0.06659613 0 0.1551559 0.05 > binom.agresti.coull(0,25) method x n mean lower upper 1 agresti-coull 0 25 0 -0.02439494 0.1575872 > ptilde(0,25,(qnorm(0.95+.05/2)))  0.06659613 Notice the mean are identical (the output from ptilde) As an aside, I do not know why binom.agresti.coull doesn't report the new mean?? I calculated as per the Wikipedia article. At least as compared to the "exact" method and N=25, binom.bayes with the Agresti-Coull shape appears better. Any thoughts on this? Question 2 follows: require(binom) ans0=c() ans1=c() for(n in (1:50)){ for(x in c(0,n)){ b=binom.exact(x,n); bayes=binom.bayes(x,n) clevel=1-pnorm(2,lower.tail = F)*2 z=2 n.tilde=n+z^2 p.tilde=1/n.tilde*(x+1/2*z^2) #see https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Agresti-Coull_Interval m=mean(c(b$lower,b$upper)); w=(m*n-x)/(0.5-m) check=(w *0.5+x)/(w+n) cat(x,'/',n,' "exact".center=',m,' bayes.mean=',bayes$mean,' ac.mean=',p.tilde,' w=',w,' weighted.avg=',check,"\n")
if(x==0) {
ans0[n]=w;
}else {
ans1[n]=w;
}
}
}
View(data.frame(xEq0=ans0,xEq1=ans1))

Output:

0 / 1  "exact".center= 0.4875  bayes.mean= 0.25  ac.mean= 0.4  w= 39  weighted.avg= 0.4875
0 / 2  "exact".center= 0.4209431  bayes.mean= 0.1666667  ac.mean= 0.3333333  w= 10.64911  weighted.avg= 0.4209431
0 / 3  "exact".center= 0.3537991  bayes.mean= 0.125  ac.mean= 0.2857143  w= 7.259856  weighted.avg= 0.3537991
0 / 4  "exact".center= 0.3011823  bayes.mean= 0.1  ac.mean= 0.25  w= 6.059467  weighted.avg= 0.3011823
0 / 5  "exact".center= 0.2609119  bayes.mean= 0.08333333  ac.mean= 0.2222222  w= 5.456396  weighted.avg= 0.2609119
0 / 6  "exact".center= 0.2296291  bayes.mean= 0.07142857  ac.mean= 0.2  w= 5.095867  weighted.avg= 0.2296291
0 / 7  "exact".center= 0.2048082  bayes.mean= 0.0625  ac.mean= 0.1818182  w= 4.856698  weighted.avg= 0.2048082
0 / 8  "exact".center= 0.1847083  bayes.mean= 0.05555556  ac.mean= 0.1666667  w= 4.686665  weighted.avg= 0.1847083
0 / 9  "exact".center= 0.1681336  bayes.mean= 0.05  ac.mean= 0.1538462  w= 4.559672  weighted.avg= 0.1681336
0 / 10  "exact".center= 0.1542486  bayes.mean= 0.04545455  ac.mean= 0.1428571  w= 4.461255  weighted.avg= 0.1542486
0 / 11  "exact".center= 0.1424571  bayes.mean= 0.04166667  ac.mean= 0.1333333  w= 4.382768  weighted.avg= 0.1424571
0 / 12  "exact".center= 0.1323242  bayes.mean= 0.03846154  ac.mean= 0.125  w= 4.318726  weighted.avg= 0.1323242
0 / 13  "exact".center= 0.1235263  bayes.mean= 0.03571429  ac.mean= 0.1176471  w= 4.265483  weighted.avg= 0.1235263
0 / 14  "exact".center= 0.1158179  bayes.mean= 0.03333333  ac.mean= 0.1111111  w= 4.220525  weighted.avg= 0.1158179

Question #2 For small values of the ideal weight parameter may be as high as 39. Would it be best to modify Agresti-Coull to adjust it's method of weight to match accordingly?

References:

Agresti-Coull Interval

Question #1: Setting $Z = 1$ with the Agresti Coull method will yield the same posterior distribution as would the Jeffreys Prior. Observe that when $Z = 1$, $Z^2 = 1$, so $\tilde n = n + 1$ and $\tilde p = \frac {X + .5} {n + 1}$.