In the continuing saga, Joseph Coveney <[email protected]> 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
[email protected]
------------------------------begin bincoverage.ado---------------------------
*! 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
--------------------------------end bincoverage.ado---------------------------
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
