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;