*  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .
*  File:  	ciprop.SPS .
*  Date:  	23-Nov-2001 .
*  Author:  Bruce Weaver, weaverb@mcmaster.ca .
*  Notes:	Get confidence interval binomial proportion using:
		- inverse Beta function
		- Ghosh (1979) method
		- Wilson score method
*  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .

* The data used here are from Table I in Newcombe (1998), Statistics
   in Medicine, Vol 17, 857-872.

* Inverse beta function method.
* This code is based on a newsgroup post by David L. Turner.
* Turner's post can be found in Rich Ulrich's Stats FAQ here:
	http://www.pitt.edu/~wpilib/statfaq/97nonpar.html .

DATA LIST LIST /x(f8.0) n(f8.0) confid(f5.3) .
BEGIN DATA.
81 263 .95
15 148 .95
0   20 .95
1   29 .95
81 263 .90
15 148 .90
0   20 .90
1   29 .90
81 263 .99
15 148 .99
0   20 .99
1   29 .99
END DATA.

compute p = x/n.
do if (x > 0).
-	compute lower1 = idf.beta((1-confid)/2,x,(n-x+1)).
else if (x = 0).
-	compute lower1 = 0.
end if.
compute upper1 = idf.beta(1-(1-confid)/2,(x+1),(n-2)).
exe.

formats p lower1 upper1 (f5.4).
list x n p lower1 upper1 .

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

* Now use method of Ghosh (1979), as described in Glass & Hopkins (1996, p 326).

* This method is described there as the "method of choice for all values
   of p and n" .

compute p = x/n.
compute q = 1-p.
compute z = probit(1-(1-confid)/2).
compute #a = (n / (n + z**2)) .
compute #b = p + (z**2/(2*n)).
compute #c = (p*q)/n.
compute #d = z**2/(4*n**2).
compute lower2 = #a * (#b - z*sqrt(#c+#d)).
compute upper2 = #a * (#b + z*sqrt(#c+#d)).
exe.

formats lower2 upper2 (f5.4).
list x n p lower2 upper2 .

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

* Wilson score method (Method 3 in Newcombe, 1998) .
* Code adapted from Robert Newcombe's code posted here:
    http://archive.uwcm.ac.uk/uwcm/ms/Robert2.html .

compute z = probit(1-(1-confid)/2).
compute p = x/n.
compute q = 1-p.
COMPUTE #x1 = 2*n*p+z**2 .
COMPUTE #x2 = z*(z**2+4*n*p*(1-p))**0.5 .
COMPUTE #x3 = 2*(n+z**2) .
COMPUTE lower3 = (#x1 - #x2) / #x3 .
COMPUTE upper3 = (#x1 + #x2) / #x3 .
EXECUTE .

formats lower3 upper3 (f5.4).
list x n p lower3 upper3.

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

* List results of all 3 methods .

list var x n confid p lower1 upper1 lower2 upper2 lower3 upper3 .

* Method 1:  Inverse beta function .
* Method 2:  Wilson score method (from Newcombe paper).
* Method 3:  Ghosh (1979) method (from Glass & Hopkins).

* Data from Newcombe (1998), Table I.
* For these data, Methods 2 and 3 produce identical results.
* Robert Newcombe looked at this, and confirmed that the Ghosh and 
   Wilson methods are algebraically identical.

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