Re: st: manually calculating confidence interval for a proportion

[Steve Samuels <sjsamuels@gmail.com>]

Correction: the problem is demonstrated with the binomial() function.
Steve

The  lower tail for -cii- seems to be in error:

. cii 484 157
. di binomial(484, 157, r(lb))


I've noted the discrepancy between the results of
. cii 484 157
and 
. di  invbinomial(484, 157, .975)
. di  invbinomial(484, 157, .025)

The  lower tail for -cii- seems to be in error:

. cii 484 157
. di invbinomial(484, 157, r(lb))

I've reported this to Stata support.

Steve




Binomial CIs of the form  Estimate +/- Bound are based on the central limit theorem  and so use  normal percentiles.

. scalar calculated_95ub_bin= 157/484 + invnormal(0.025)*(sqrt((1/484)*(157/484)*(1-(157/484))))
. scalar calculated_95lb_bin= 157/484 - invnormal(0.025)*(sqrt((1/484)*(157/484)*(1-(157/484))))

** Note use of 0.025

. dis calculated_95ub_bin
. dis calculated_95lb_bin
. cii 484 157

** The exact CIs are calculated as follows:

. di  invbinomial(484, 157, .975)
. di  invbinomial(484, 157, .025)

In the absence of the invbinomial() function, the exact CIs can be calculated in
terms of the F distribution or the beta distribution.  See the Wikipedia
entry http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval


Steve
sjsamuels@gmail.com

On Jul 12, 2013, at 2:49 PM, Wood, Mollie Elizabeth wrote:

I am putting together a teaching example for a beginning programming lab, in which I ask students to manually calculate the upper and lower confidence limits for a variety of point estimates, and compare their results to the stata commands (ci, cii).  However, I haven't had much look making it work myself!
> 
I used this code to obtain 90% CIs for a mean:
scalar calculated_90ub = r(mean) + invttail(_N - 1,0.05)*LOS_se
scalar calculated_90lb = r(mean) - invttail(_N - 1,0.05)*LOS_se
dis calculated_90ub
dis calculated_90lb
> 
Which worked well, and matched the output of Stata's command, ci [var], level(90)
> 
However, when I tried the below command for the 95% CI around a proportion, the results are quite different the stata command, cii N K, level(95)
> 
scalar calculated_95ub_bin= 157/484 + invbinomial(484,157,.05)*(sqrt((1/484)*(157/484)*(1-(157/484))))
scalar calculated_95lb_bin= 157/484 - invbinomial(484,157,.05)*(sqrt((1/484)*(157/484)*(1-(157/484))))
> 
dis calculated_95ub_bin
dis calculated_95lb_bin
> 
I did experiment with several other related functions (invbinomialtail, etc), none of which accurately reproduced the output from cii.
> 
If anyone has a suggestion as to how I'm going wrong, I'd appreciate it very much!


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/