proc iml;

 /*-------------------------------------------------------------*/
 /*                 The MULTINOM Module                         */
 /*                                                             */
 /* This module constructs large-sample confidence limits for   */
 /* various linear combinations of the K multinomial            */
 /* proportions.                                                */
 /*                                                             */
 /* Input:                                                      */
 /*    1) NVEC, a 1xK vector containing the cell sizes.         */
 /*    2) TYPE, a character string that specifies whether you   */
 /*       want the confidence limits to control individual or   */
 /*       overall error rate.  The choices are "i"              */
 /*       (for individual) or "s" (for simultaneous).           */
 /*    3) ALPHA, the type I error rate.                         */
 /*    4) CONT, a MxK matrix containing the M linear            */
 /*       combinations of interest.  If CONT is set to missing, */
 /*       MULTINOM automatically constructs confidence          */
 /*       intervals for all K(K-1)/2 pairwise differences       */
 /*-------------------------------------------------------------*/

start multinom(nvec,type,alpha,cont);
reset noname fw=2;
type=upcase(type);
if min(nvec)<0 | sum(abs(round(nvec,1)-nvec))>1e-4 then do;
   print "ERROR: Cell Sizes Must Be Nonnegative Integers";
   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);
if cont^=. & ncol(cont)^=k then do;
   print "ERROR: Number of Columns in CONT Must Equal K";
   return;
end;
n=sum(nvec);
pvec=nvec/n;
varcov=(diag(pvec)-pvec`*pvec)/n;
cl=(1-alpha)#100;
c2="Estimate" "LCL" "UCL";
if type='I' then label='Individual';
if type='S' then label='Simultaneous';
if cont=. then do;
   c1="I" "J";
   print "All Pairwise Differences";
   print cl "%" label "Large-Sample Confidence Limits For Pi-Pj";
   c=round(k#(k-1)/2,1);
   contr=j(c,k,0);
   i=1;  j=2;
   do row = 1 to c;
      index=index//(i||j);
      contr[row,i]=1;
      contr[row,j]=-1;
      j=j+1;
      if j>k then do;
         i=i+1;
         j=i+1;
      end;
   end;
end;
else do;
   c1="Ci";
   print "User-Defined Linear Combinations";
   print cl "%" label "Large-Sample Confidence Limits";
   contr=cont;
   c=nrow(contr);
   index=cont;
end;
if type='I' then z=probit(1-alpha/2);
if type='S' then z=probit(1-alpha/(2#c));
do i = 1 to c;
   ci=contr[i,];
   mean=ci*pvec`;
   var=ci*varcov*ci`;
   inc=z#sqrt(var);
   lcl=mean-inc;
   ucl=mean+inc;
   table=table//(mean||lcl||ucl);
end;
print index[colname=c1 format=2.0] table[colname=c2 format=6.4];
finish multinom;


 /*-------------------------------------------------------------*/
 /*                 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;