1
$\begingroup$

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?

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)))
[1] 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

$\endgroup$
1
$\begingroup$

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}$.

Question #2: In general, the prior should not change based on the sample size. The prior should reflect all the information that is known before any information has been observed. Whether a small or large amount of data will be observed, the prior knowledge remains the same.

$\endgroup$
  • $\begingroup$ In regards to Question #1, the question is what are the Beta distribution shape parameters to match Agresti Coull rather than what Agresti Coull parameters to match Jeffrey's Prior. I still awarded an up vote however because what you wrote is insightful. $\endgroup$ – Chris Aug 2 '16 at 1:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.