  # Re: st: Binomial confidence intervals

 From rgutierrez@stata.com (Roberto G. Gutierrez, StataCorp) To statalist@hsphsun2.harvard.edu Subject Re: st: Binomial confidence intervals Date Wed, 08 Sep 2004 14:21:40 -0500

```In the continuing saga, Joseph Coveney <jcoveney@bigplanet.com> summarizes
the results of his simulation study.

> In the absence of access to the article, you can run -simulate- calling a
> program such as the -exbinci- ditty in the do-file below, and make a choice
> suitable to your circumstances based on the results.  I wrote the do-file
> below in an attempt to illustrate Bobby Gutierrez's point to the list.  In
> order to run it, you'll need to install Joseph Hilbe's -rnd- suite from SSC.

> In the do-file below, with 10 trials and a population mean of 50% (these are
> options in the program that you can change to suit your circumstances), the
> true parameter lies within the 95% confidence interval 9797 times out of
> 10000 experiments for each of the methods.  This compares with a 95%
> confidence interval's expectation to contain the parameter 9500 times out of
> the 10000 experiments.  (A 95% confidence interval is supposed to contain
> the population parameter 95% of the time over the long run.)

[...]

As an alternative to simulation, Nick Cox and I offer the following code for
-bincoverage-, which will calculate the coverage probabilities exactly (and
here I do mean "exactly").

. bincoverage, n(400) p(0.21)

For N = 400 and p = .21, the true coverage probability of the nominal
95% exact CI is 0.9574.

. bincoverage, n(400) p(0.21) jeffreys

For N = 400 and p = .21, the true coverage probability of the nominal
95% jeffreys CI is 0.9505.

etc.

--Bobby
rgutierrez@stata.com

*! version 1.0.4  08sep2004  rgg/njc
program bincoverage, rclass
version 8.2
syntax [, n(int 10) p(real 0.5) Level(int `c(level)') ///
EXAct Wilson Agresti Jeffreys]

local flavor `exact' `wilson' `agresti' `jeffreys'
local m : word count `flavor'

if `m' > 1 {
di as err "{p}You may only specify one of exact, wilson, " _c
di as err "agresti, or jeffreys{p_end}"
exit 198
}
else if `m' == 0 {
local flavor exact
}
local mid = round(`n' * `p',1)
local kmin `mid'
local kmax `mid'

// begin in the middle, branch out until you get no coverage

while `kmin' >= 0 {
qui cii `n' `kmin', `flavor' level(`level')
if missing(r(lb)) | missing(r(ub)) {
local miss miss
}
if r(lb) < `p' & r(ub) > `p' {
local --kmin
}
else {
continue, break
}
}
while `kmax' <= `n' {
qui cii `n' `kmax', `flavor' level(`level')
if missing(r(lb)) | missing(r(ub)) {
local miss miss
}
if r(lb) < `p' & r(ub) > `p' {
local ++kmax
}
else {
continue, break
}
}
if "`miss'" != "" {
di as err "confidence limits were missing in one " _c
di as err "or more outcomes "
di as err "n is probably too large"
exit 498
}
tempname cp
scalar `cp' = cond(`kmin'==`kmax',0, ///
Binomial(`n',`kmin'+1,`p') - Binomial(`n',`kmax',`p'))
di as txt _n "For N = " as res `n'  as txt " and p = " as res `p' ///
as txt ", the true coverage probability of the nominal" ///
_n "    " ///
as res "`level'% `flavor'" ///
as txt " CI is " as res %6.4f `cp' ///
as txt "."

return scalar n = `n'
return scalar p = `p'
return local type `flavor'
return scalar coverage = `cp'
end