options ps=65 ls=65; *This is power.pilot.sas; *Two data files are read by this program. The first is read into rawdata and contains a group indicator variable and scores on p variables. The second is read into param and contains the group size to be considered for the study, alpha, and the number of replications for the simulation.; data rawdata; input group var1-var3; datalines; 1 1 6 3 1 2 7 3 1 3 7 1 1 9 4 4 1 6 6 5 1 6 4 3 1 3 9 6 1 4 9 7 1 5 4 3 1 4 10 8 2 1 8 4 2 0 6 0 2 2 8 3 2 2 3 0 2 6 7 5 2 1 10 4 2 3 9 3 2 1 11 6 2 5 8 3 2 0 7 2 proc sort; by group; proc corr cov data=rawdata outp=cov; by group; var var1-var3; data param; input n alpha nreps; datalines; 60 .05 100 proc iml; edit cov; read all into temp; p=ncol(temp)-1; cov1=temp[1:p,2:p+1]; cov2=temp[2*p+4:3*p+3,2:p+1]; n1=temp[p+3,2]; n2=temp[3*p+6,2]; cov=((n1-1)*cov1+(n2-1)*cov2)/(n1 + n2 - 2); barx1=temp[p+1,2:p+1]; barx2=temp[3*p+4,2:p+1]; d=(barx1-barx2)*(inv(sqrt(diag(cov)))); print barx1, barx2, cov1, cov2, cov,d; edit param; read all into temp1; n=temp1[1,1]; df=2*(n-1); alpha=temp1[1,2]; nreps=temp1[1,3]; 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; bx1=barx1+((rannor(tbx))*U)/(sqrt(n)); bx2=barx2+((rannor(tbx))*U)/(sqrt(n)); t=abs(((sqrt(n/2))*(bx1-bx2)*(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;