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" "   ChiSq"};
if type='I' then label='Individual';
if type='S' then label='Simultaneous';
if cont=. then do;
   c1={"I" "J"};
   title1="All Pairwise Differences";
   title2=concat(char(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";
   title1="User-Defined Linear Combinations";
   title2=concat(char(cl),"% ",label," Large-Sample Confidence Limits");
   contr=cont;
   c=nrow(contr);
   index=rowcat(char(cont));
end;
if type='I' then r=1; else r=c;
do i = 1 to c;
   ci=contr[i,];
   mean=ci*pvec`;
   var=ci*varcov*ci`;
   z=probit(1-alpha/(2#r));
   inc=z#sqrt(var);
   lcl=mean-inc;
   ucl=mean+inc;
   if mean ^= 0 & var ^= 0 then do;
     wald=mean*mean/var; 
     pval=pval//(r#(1-probchi(wald,1)));
   end;
   else do;
     wald=0; pval=1;
   end;
   table=table//(mean||lcl||ucl||wald);
end;
print /;
title "The MULTINOM module";
print title1, title2;
print index[colname=c1] 
      table[colname=c2 format=6.4]
       pval[colname="Pr>ChiSq" format=pvalue.];
finish multinom;