[R] Confidence interval for binomial variance

[Robert A. LaBudde]

Based on simulations, I've come up with a simple function to compute 
the confidence interval for the variance of the binomial variance, 
where the true variance is

	v	= rho*(1-rho)/n

where rho = true probability of success and n = # of trials.

For x = # successes observed in n trials, p = x / n as usual.

For p < 0.25 or p > 0.75, I use the proportion-based transformed 
confidence interval, as suggested by Moshe Olshansky. I used the 
Wilson interval, but some might prefer the Blaker interval. For p >= 
0.25 or p <= 0.75, I use the standard chi-square based confidence 
interval (which is very conservative for this problem in this range). 
The composite method gives reliable results over a wide range of rho.

This should suit the purpose until someone comes up with a more 
theoretically sound method. Bootstrap is not reliable for n < 50, at least.

#09.26.08 02.50 binomVarCI.r
#copyright 2008 by Robert A LaBudde, all rights reserved
#CI for binomial sample variance
#created: 09.26.08 by r.a. labudde
#changes:

require('binGroup')

binomVarCI<- function (n, x, conf=0.95) {
   p<- x/n #proportion
   if (p<0.25 | p>0.75 | x==0 | x==n) {  #use proportion-based CI
     pCI<- binWilson(n, x, conf.level=conf)  #CI for proportion
     vCI<- sort(c(pCI[1]*(1-pCI[1])/(n-1), pCI[2]*(1-pCI[2])/(n-1)))
   } else {  #use chi-square-based CI
     phiL<- qchisq(0.025, n-1)/(n-1)
     phiU<- qchisq(0.975, n-1)/(n-1)
     #vest<- p*(1-p)/(n-1))  #variance estimate
     vCI<- c(vest/phiU, vest/phiL) #chi-square-based
   }
   return (vCI)
}

Here is a test program to measure coverage:

#09.26.08 03.10 tbinomVarCI.r
#copyright 2008 by Robert A LaBudde, all rights reserved
#test of CI for binomial sample variance
#created: 09.26.08 by r.a. labudde
#changes:

nReal <- 1000

for (POD in c(0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.50)) {
   cat('\nPOD: ', sprintf('%8.4f', POD), '\n')
   for (nRepl in c(6, 12, 20, 50)) {
     vtrue<- POD*(1-POD)/nRepl
     pcover<- 0
     for (iReal in 1:nReal) {
       x<- rbinom(1, nRepl, POD)
       vCI<- binomVarCI(nRepl, x)
       #vest<- (x/nRepl)*(1 - x/nRepl)/(nRepl-1)
       if (vtrue >= vCI[1] & vtrue<= vCI[2]) pcover<- pcover + 1
     } #iReal
     pcover<- pcover/nReal
     cat('n: ', sprintf('%4i', nRepl), ' Coverage: ', 
sprintf('%8.4f', pcover), '\n')
   } #nRepl
} #POD
================================================================
Robert A. LaBudde, PhD, PAS, Dpl. ACAFS  e-mail: ral at lcfltd.com
Least Cost Formulations, Ltd.            URL: http://lcfltd.com/
824 Timberlake Drive                     Tel: 757-467-0954
Virginia Beach, VA 23464-3239            Fax: 757-467-2947

"Vere scire est per causas scire"