Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

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

------------------------------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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index