calculating true confidence intervals across binomial distributions

[Jack McDonnell]

I'm attempting to reproduce in Stata a technique from a statistical coding textbook ("Analysis of Categorical Data with R", Bilder and Loughlin).

I ultimately want to show visually how near confidence intervals for a binomial distribution (using the Wald technique) approach the idealized 95% for a range of probability values - 0.001(0.0005)0.999). In other words, something like the graphic attached below (which was generated in R using the author's code).

I'm really trying to iterate the sample code below across the range of 0.001(0.0005)0.999, whereas the below code only gets me pi=0.001 and I would like to get it for the whole range.

clear
set obs 1998
gen w = rbinomial(40,0.001)
gen pi_hat = w/40
gen var_wald = (pi_hat*(1-pi_hat)/40)
gen lower = pi_hat + (invnormal(.025) * sqrt(var_wald))
gen upper = pi_hat - (invnormal(.025) * sqrt(var_wald))
gen pmf=0
replace pmf=binomialp(40,w,pi)
gen save = 1 if pi>lower & pi<upper

gen CI = save/1998

The iterative techniques I've tried encompass many variations of the following, but I invariably get endless loops or broken code that doesn't give me what I want.

forvalues j = 1/1998{
forvalues i = 0.001(0.0005)0.999{
replace pi = `i'
replace w = rbinomial(40,pi)
replace pi_hat = w/40
replace var_wald = (pi_hat*(1-pi_hat)/40)
replace lower = pi_hat + (invnormal(.025) * sqrt(var_wald))
replace upper = pi_hat - (invnormal(.025) * sqrt(var_wald))
replace pmf=0
replace pmf=binomialp(40,w,pi)
replace save = 0
replace save = 1 if pi>lower & pi<upper
}
replace CI = save/1998 in `j'
}

My goal is for the variable CI to contain a list of all the unique proportions describing how often the Wald confidence interval contained the "true" pi, for a given value of pi.

This is my first time posting in the group, although I've followed for years to learn the program better. I do apologize if I've made mistakes in formulating this question.

Thanks in advance,


Jack McDonnell
[ATTACH=CONFIG]n1555902[/ATTACH]

[Jack McDonnell]

the graphic I referenced

[Mike Lacy]

The way I've done this sort of thing is to write a program that simulates the CI for a given N, PI, etc., and then use it within the Stata suite of -postfile- commands to accumulate the various results into a data set. When that is done, you can graph etc. with the results. See -help postfile-, which has useful examples. Here's an example program and a sketch of what you'd do with -postfile-.

// Quick and dirty program to simulate the CI performance of Wald technique
// for given PI, samplesize, confidence level, and number of simulation
// repetitions.  It returns the coverage results in the r() return list.
cap prog drop checkCI
prog checkCI, rclass
syntax, pi(real) n(real) level(real) reps(real)
clear
quiet set obs `n'
tempvar y
quiet gen byte `y' = .
local covers = 0
forval i = 1/`reps' {
   // n Bernouilli trials
   quiet replace `y' = runiform() < `pi'
   // pi_hat and CI from this rep
   quiet ci proportions `y', wald level(`level')  // let Stata do the CI, not you
   if inrange(`pi', r(lb), r(ub)) {
      local ++covers
   }
}
return scalar nreps = `reps'
return scalar ncovers = `covers'
return scalar pcovers = `covers'/`reps'
end
//
// An example use of the preceding program
checkCI, pi(0.5) level(95) reps(1000) n(30)
return list

With this program, you can set up to loop over your PI values, and accumulate the results. In outline form, you'd have something like:

postfile .....using YourResultsFile.dta   // set up postfile
set seed 75654
forvalues p = ..... your choice of PI values
   checkCI, pi(`p') level .....
   post .... // post results for this particular pi value
}
postclose .....
clear
use YourResultsFile.dta

[Jack McDonnell]

thanks Mike, this is really helpful