options ps=65 ls=65;
*This is power.statistic.sas;
*Three data files are read into this program.  The first is
read into rs and contains a correlation matrix for p variables.
The second is read into ds and contains effect sizes on the
same three variables.  The third is read into param and
contains the sample size to be considered for the study,
alpha, and the number of replications.;
data rs;
input var1-var3;
datalines;
1.00  .25  .50
 .25 1.00  .33
 .50  .33 1.00
data ds;
input d1-d3;
datalines;
.5 .5 .2
data param;
input n alpha nreps;
datalines;
60 .05 100
proc iml;
edit rs;
read all into corr;
print corr;
edit ds;
read all into d;
print d;
edit param;
read all into temp1;
n=temp1[1,1];
alpha=temp1[1,2];
nreps = temp1[1,3];
p=nrow(corr);
ptest=ncol(d);
if p ^= ptest then do;
print '                        ERROR
';
print 'number of effect sizes and number of variables are not
equal';
end;
do ii = 1 to p;
do jj = ii+1 to p;
if corr[ii,jj] ^= corr[jj,ii] then do;
print '                        ERROR
';
print 'the  correlation matrix is not symmetric';
end;
end;
end;
barxd=((sqrt(diag(corr)))*d`)`;
barxd=barxd`;
df=2*(n-1);
nt=p;
sigibh=j(nt,1,0);
sigi  =j(nt,1,0);
sigall=0;
sigallbh=0;
siganybh=0;
sigany  =0;
U = ROOT(Corr);
tbx=j(1,p,0);
do reps = 1 to nreps;
   ts1=j(p,p,0);
     do ii=1 to p-1;
     do jj=ii+1 to p;
        ts1[jj,ii]=rannor(0);
     end;
     end;
     do ii=1 to p;
        ts1[ii,ii]=sqrt(2*rangam(0,(n-ii)/2));
     end;
   stemp1=(ts1*ts1`);
   ts2=j(p,p,0);
     do ii=1 to p-1;
     do jj=ii+1 to p;
        ts2[jj,ii]=rannor(0);
     end;
     end;
     do ii=1 to p;
        ts2[ii,ii]=sqrt(2*rangam(0,(n-ii)/2));
     end;
   stemp2=(ts2*ts2`);
   stemp=(stemp1+stemp2)/(2*n-2);
   stemp=u`*stemp*u;*following for debugging;
   bxd=barxd`+((rannor(tbx))*U)*(sqrt(2/n));
   t=abs(((sqrt(n/2))*(bxd)*(inv(sqrt(diag(stemp))))))`;
   rankt=rank(t);
   allcnt=0;
   comp=0;
   do mm=1 to nt;
   do nn=1 to nt;
      if rankt[nn,1] = mm then do;
         if (t[nn,1]) > tinv(1-alpha/(2*mm),2*n-2) then do;
            comp=t[nn,1];
            sigany  =sigany  +1;
            if mm=1 then allcnt=allcnt+p;
            do rr=1 to nt;
               if t[rr,1] >= comp then
               sigi[rr,1]=sigi[rr,1]+1;
            end;
            nn=nt;
            mm=nt;
         end;
      end;
   end;
   end;
   If allcnt=p then sigall=sigall+1;
   anycount=0;
   allcntbh=0;
   do oo=1 to nt;
   do qq=1 to nt;
      if rankt[qq,1]=nt-oo+1 then do;
         if (t[qq,1]) > tinv(1-alpha/(2*(nt-oo+1)),2*n-2) then
do;
            sigibh[qq,1]=sigibh[qq,1]+1;
            anycount=anycount+1;
            allcntbh=allcntbh+1;
         end;
         if (t[qq,1]) <= tinv(1-alpha/(2*(nt-oo+1)),2*n-2) then
do;
            qq=nt;
            oo=nt;
         end;
      end;
   end;
   end;
if anycount > 0 then siganybh=siganybh+1;
If allcntbh=p then sigallbh=sigallbh+1;
end;
Print 'The number of replications is';
print nreps;
Print 'The group size is';
print n;
Print 'Note.  N is the number of Ss in each group.';
Print 'Following results for Hochberg';
anypower  =sigany  /(reps-1);
avgpower  =sigi  [+]/((reps-1)*nt);
allpower  =sigall  /(reps-1);
indpower=sigi/(reps-1);
print anypower, avgpower, allpower, indpower;
Print 'Following results for Bonferroni-Holm';
anypower=siganybh/(reps-1);
avgpower=sigibh[+]/((reps-1)*nt);
allpower=sigallbh/(reps-1);
indpower=sigibh/(reps-1);
print anypower, avgpower, allpower, indpower;
quit;