/* Appendix 2                                              */

/* Functionally, these modified IML modules are the same   */
/* as the original modules(*) except for the addition of   */
/* the randomization test to test for differences between  */
/* to test for differences between the means of two        */
/* independent samples. An example that illustrates the    */
/* use of the PERMTEST module for this addition is shown   */
/* at the end of the Appendix.                             */
/*                                                         */
/* (*) See "A PROC IML Program to obtain Exact             */
/* Significance Levels in the Nonparametric                */
/* Two-Independent-Samples: Location Problem" in           */
/* Observations, Fourth Quarter 1993.                      */

proc iml;

   /*          The Choose Module                           */
   /* We consider all the possible choices of M objects    */
   /* from a set of M+N. There is a sequential way of      */
   /* ordering them. The module CHOOSE gives you the next  */
   /* choice when you supply a given selection and the     */
   /* parameters M and N.  The vector VEC is input and     */
   /* consists of integers from 1 to M+N. Examples if you  */
   /* are chosing five objects from a set of 8:            */
   /* 1) Input mvec0=[1 2 3 5 7] and get mvec1=[1 2 3 5 8] */
   /* 2) Input mvec0=[1 3 4 7 8] and get mvec1=[1 3 5 6 7] */
   /* This module is used inside the module PERMTEST.      */

start choose(vec,mvec0,m,n, mvec1);
mvec1=mvec0;
do place=m to 1 by -1;
   if mvec0[1,place]<n+place then do;
      mvec1[1,place]=mvec1[1,place]+1;
      if place<m then mvec1[1,place+1:m]=
                 vec[1,place+1:m]+mvec1[1,place]-place;
      return;
   end;
end;
finish;

/*             The PERMTEST Module                         */
/*                                                         */
/* This module performs the permutation test.              */
/*                                                         */
/* Input:                                                  */
/*  1) TEST, a character string that specifies the         */ 
/*     scores used.  The possible choices are MEAN,        */ 
/*     WILCOXON, MEDIAN, VW, and SAVAGE.                   */ 
/*  2) M, the sample size of the CLASS variable level      */ 
/*     designated "first".                                 */ 
/*  3) N, the sample size of the CLASS variable level      */ 
/*     designated "second".                                */ 
/*  4) SUMSCORE, the sum of the scores for the "first"     */ 
/*     class level for the desired test.                   */ 
/*  5) YVEC, a (M+N) x 1 vector containing the             */ 
/*     dependent variable values. Only used when the       */ 
/*     choice for TEST is "MEAN".                          */ 
/*                                                         */ 
/* Returns PVAL, a 1x3 vector that gives the exact         */ 
/* significance level for lower-tailed, upper-tailed,      */ 
/* and two-tailed hypotheses.                              */ 

start permtest(test,m,n,sumscore,yvec, pval);
test=upcase(test);
mn=m+n;
size=round(gamma(mn+1)/(gamma(m+1)#gamma(n+1)),1);
vec=1:mn;
pval=j(1,3,0);
mvec0=vec[1,1:m];
if test='MEAN' then ysum=yvec[+];
if test='SAVAGE' then savvec=1/(mn-vec+1);
do i = 1 to size;
   if test='.' | test='WILCOXON' then score0=mvec0;
   if test='MEDIAN' then score0=(mvec0>(mn+1)/2);
   if test='VW' then score0=probit(mvec0/(mn+1));
   if test='SAVAGE' then do;
      score0=j(1,m,.);
      do j = 1 to m;
         mj=mvec0[1,j];
         score0[1,j]=sum(savvec[1,1:mj])-1;
       end;
   end;
   if test='MEAN' then score0=yvec[mvec0,];
   rankstat=sum(score0);
   if rankstat<=sumscore then pval[1,1]=pval[1,1]+1;
   if rankstat>=sumscore then pval[1,2]=pval[1,2]+1;
   if test='MEAN' then do;
      if abs(mn#rankstat-m#ysum)>=abs(mn#sumscore-m#ysum)
      then pval[1,3]=pval[1,3]+1;
   end;
   if i<size then do;
      run choose(vec,mvec0,m,n,mvec1);
      mvec0=mvec1;
   end;
end;
if test^='MEAN' then pval[1,3]=2*min(pval[1,1],pval[1,2]);
pval=pval/size;
if pval[1,3]>1 then pval[1,3]=1;
reset noname;
chead={"Lower-Tail" "Upper-Tail" "Two-Tail"};
print pval[colname=chead];
finish;


   /* For Example 2                                        */
   /* Use 'Mean' for Test                                  */
   /* M=4 N=6                                              */
   /* SUMSCORE = Sum of Y's from sample 1 = 50.8           */
yvec={12.6 14.4 12.3 11.4 8.3 12.6 8.4 11.5 10.3 6.2}`;
run permtest('mean',4,6,50.8,yvec, pval);