*  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .
*  File:  	clopper-pearson.SPS .
*  Date:  	23-Nov-2001 .
*  Author:  Bruce Weaver, weaverb@mcmaster.ca .
*  Notes:	Get confidence interval binomial proportion using the
		 Clopper-Pearson "exact" method .
*  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .

* X = the observed number of events.
* N = the total number of trials.
* Enter X and N on the two lines flagged below.
* Enter the confidence level on the following line (e.g., 0.95) .

new file.
input program.
loop #i = 1 to 10000.
	compute x = 81.		/*  Enter X on this line  */
	compute n = 263.		/*  Enter N on this line  */
	compute confid = 0.95.	/*  Enter confidence level here */
	compute pobs = x/n.
	compute alpha = 1 - confid.
	compute p = #i/10000.
	end case.
end loop.
end file.
end input program.
exe.

COMPUTE sum0toX = CDF.BINOM(x,n,p) .
compute sumXtoN = 1 - CDF.BINOM(x-1,n,p) .
exe.

* Create flags for:
	- sum of p from 0 to X < alpha/2 
	- sum of p from X to N < alpha/2 .

compute flag1 = (sum0toX < alpha/2).
compute flag2 = (sumXtoN < alpha/2).
exe.
formats flag1 flag2 (f2.0).


* Sort cases so largest FLAGGED value of sum0toX is on 1st row .
* Value of p on this row is one limit of the confidence interval for x/N .
* Flag this record as one of the CI limits.

sort cases by flag1(d) sum0toX (d).
compute flag = ($casenum=1).
exe.

* Now sort cases so largest FLAGGED value of sumXtoN is on 1st row .
* Then keep 1st row only, AND the row with FLAG=1.
* Values of p on these two rows are the limits of the 
   Clopper-Pearson "exact" confidence interval for x/N .

sort cases by flag2(d) sumXtoN (d).
if not(flag) flag = ($casenum=1).
select if flag.
exe.

* Now compute limits using the Wilson/Ghosh (1979) method .

compute q = 1-pobs.
compute z = probit(1-(1-confid)/2).
compute #a = (n / (n + z**2)) .
compute #b = pobs + (z**2/(2*n)).
compute #c = (pobs*q)/n.
compute #d = z**2/(4*n**2).
compute wglower = #a * (#b - z*sqrt(#c+#d)).
compute wgupper = #a * (#b + z*sqrt(#c+#d)).
exe.

format x n (f5.0) / alpha confid pobs p sum0toX sumXtoN wglower wgupper (f8.4).

rename var (p = limits).
var lab
 pobs 'Obsvered proportion'
 limits 'Clopper-Pearson "exact" limits'
 confid 'Confidence level'
 wglower 'Lower limit (W/G)'
 wgupper 'Upper limit (W/G)'.

summarize 
 table = x n pobs limits confid
 /title = 'Clopper-Pearson "exact" confidence interval for a binomial proportion'
 /format = nocasenum list
 /cells = none.

compute f = ($casenum=1).
filter by f.
exe.

summarize 
 table = x n pobs wglower wgupper confid
 /title = 'Confidence interval for a binomial proportion: Wilson/Ghosh method'
 /format = nocasenum list
 /cells = none.

*  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .