proc iml;

 /*-------------------------------------------------------------*/
 /*                 The MULTIBIN Module                         */
 /*                                                             */
 /* This module constructs large-sample conservative &          */
 /* optimistic confidence intervals for the pairwise            */
 /* differences between K binomial proportions when individual  */
 /* cell sizes are unknown.  Based on the two intervals, the    */
 /* the proportion differences are labeled as "significant",    */
 /* "non-significant", or "inconclusive".                       */
 /*                                                             */
 /* Input:                                                      */
 /*    1) N, the common sample size.                            */
 /*    2) NVEC, a 1xK vector containing the number of           */
 /*       "successes" for each binomial response                */
 /*    3) TYPE, a character string that specifies if you want   */
 /*       individual or simultaneous confidence inferences      */
 /*       Choices are "i" or "s".  Only the "conservative       */
 /*       intervals are adjusted to be simultaneous.            */
 /*    4) ALPHA, the type I error rate.                         */
 /*-------------------------------------------------------------*/

start multibin(n,nvec,type,alpha);
reset noname;
type=upcase(type);
if n<=0 | abs(round(n,1)-n)>1e-4 then do;
   print "ERROR: N Must Be A Positive Noninteger";
   return;
end;
if min(nvec)<0 | max(nvec)>n |
   sum(abs(round(nvec,1)-nvec))>1e-4 then do;
   print "ERROR: NVEC Values Must be Integers & in [0,N]";
   return;
end;
if type^="I" & type^="S" then do;
   print "ERROR: TYPE Parameter Must Be 'I' or 'S'";
   return;
end;
if alpha<0 | alpha>1 then do;
   print "ERROR: ALPHA Must be in [0,1]";
   return;
end;
k=ncol(nvec);
pvec=nvec/n;
qvec=1-pvec;
c=round((k#(k-1)/2),1);
cl=100#(1-alpha);
zo=probit(1-alpha/2);
if type='I' then do;
   zc=zo;
   label='Individual';
end;
if type='S' then do;
   zc=probit(1-alpha/(2#c));
   label='Simultaneous';
end;
do i = 1 to k-1;
   pi=pvec[1,i];
   qi=qvec[1,i];
   do j = i+1 to k;
      result="  ?";
      pj=pvec[1,j];
      qj=qvec[1,j];
      mean=pi-pj;
      cvar=(1/n)#(min(pi+pj,qi+qj)-(pi-pj)##2);
      cinc=zc#sqrt(cvar);
      clcl=mean-cinc;
      cucl=mean+cinc;
      if clcl>0 | cucl<0 then result="  S";
      diff=abs(pi-pj);
      ovar=(1/n)#diff#(1-diff);
      oinc=zo#sqrt(ovar);
      olcl=mean-oinc;
      oucl=mean+oinc;
      if olcl<=0 & oucl>=0 then result="  N";
      index=index//(i||j);
      table=table//(mean||clcl||cucl||olcl||oucl);
      rvec=rvec//result;
   end;
end;
ch1={"I" "J"};
ch2={"Estimate" "C LCL" "C UCL" "O LCL" "O UCL"};
ch3="Result";
print cl "%" label "Large-Sample Confidence Limits For Pi-Pj";
print index[colname=ch1 format=2.0]
      table[colname=ch2 format=6.4] rvec[colname=ch3];
finish multibin;