options ps=65 ls=65;
*This is power.pilot.statistics.sas;
*Three data files are read into this program.  The first is
read into rawdata and contains scores on 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 minimum sample size to be considered for the
study, the maximum sample size to be considered for the
study, alpha, and the target power.;
data rawdata;
input var1-var3;
datalines;
 1  6 3
 2  7 3
 3  7 1
 9  4 4
 6  6 5
 6  4 3
 3  9 6
 4  9 7
 5  4 3
 4 10 8
 1  8 4
 0  6 0
 2  8 3
 2  3 0
 6  7 5
 1 10 4
 3  9 3
 1 11 6
 5  8 3
 0  7 2
proc corr data=rawdata cov outp=rs noprint;
  var var1-var3;
*following for debugging;
*proc print;
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 temp;
p=ncol(temp);
cov=temp[1:p,1:p];
print cov;
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];
ptest=ncol(d);
if p ^= ptest then do;
print '                       ERROR
'
print 'number of effect sizes and number of variables are not
equal';
minn=maxn + 1;
end;barxd=((sqrt(diag(cov)))*d`)`;
print barxd;
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(Cov);
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;
   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;