[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

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/

- Prev by Date:
**Re: st: Scheme error message** - Next by Date:
**RE: st: help: mlogit and outreg** - Previous by thread:
**RE: st: Binomial confidence intervals** - Next by thread:
**st: RE: Binomial confidence intervals** - Index(es):

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