options pagesize=65 linesize=65; data param; input alpha tpower rho maxn j k l int collapse; datalines; .05 .70 .00 200 2 2 4 2 1 data means; input mu1 mu2 mu3 mu4 mu5 mu6 mu7 mu8 mu9 mu10 mu11 mu12 mu13 mu14 mu15 mu16; datalines; .0 .2 .4 .5 .1 .6 .2 .4 .1 .3 .5 .6 .2 .7 .3 .5 proc iml; edit param; read all into x; alpha=x[1,1]; tpower=x[1,2]; rho=x[1,3]; maxn=x[1,4]; j=x[1,5]; k=x[1,6]; l=x[1,7]; test=ncol(x); if l=0 then test=6; *print test; if test=6 then do; print alpha tpower rho j k; edit means; read all into mur; mu=j(j,k,0); do jj=1 to j; do kk=1 to k; mu[jj,kk]=mur[1,kk+(jj-1)*k]; end; end; print mu; cj=i(j-1)//j(1,j-1,-1); ck=i(k-1)//j(1,k-1,-1); *print cj, ck; c=cj@ck; *print c; do n=2 to maxn; df2=j*k*(n-1); if rho^=0 then df2=df2-1; ni=inv(n*i(j*k)); lambda=(mur*c*(inv(c`*ni*c))*c`*mur`)/(1-rho**2); cval=finv(1-alpha,(j-1)*(k-1),df2); power=1-probf(cval,(j-1)*(k-1),df2,lambda); if power < tpower & n = maxn then do; comment = 'raise maxn '; power=round(power,.01); print comment n power; end; if power >= tpower then do; omegas=(1-rho**2)*lambda/(lambda+n*j*k); print omegas; cell_n=n; comment ='analysis ok'; power=round(power,.01); omegas=round(omegas,.001); print 'Power for the two-factor interaction'; print comment omegas cell_n power; n=maxn; end; end; end; if test=9 then do; l=x[1,7]; int=x[1,8]; collapse=x[1,9]; print alpha tpower rho j k l; edit means; read all into mur; *print mur; mu=j(j*k,l,0); do jj=1 to j; do kk=1 to k; do ll=1 to l; mm=kk+(jj-1)*k; mu[mm,ll]=mur[1,ll+(mm-1)*l]; end; end; end; print mu; cj=i(j-1)//j(1,j-1,-1); ck=i(k-1)//j(1,k-1,-1); cl=i(l-1)//j(1,l-1,-1); if int=2 & collapse=1 then cj=j(j,1,1); if int=2 & collapse=2 then ck=j(k,1,1); if int=2 & collapse=3 then cl=j(l,1,1); *print cj, ck, cl; c=cj@ck@cl; *print c; do n=2 to maxn; df2=j*k*l*(n-1); if rho^=0 then df2=df2-1; ni=inv(n*i(j*k*l)); lambda=(mur*c*(inv(c`*ni*c))*c`*mur`)/(1-rho**2);; cval=finv(1-alpha,(j-1)*(k-1)*(l-1),df2); power=1-probf(cval,(j-1)*(k-1)*(l-1),df2,lambda); if power < tpower & n = maxn then do; comment = 'raise maxn '; power=round(power,.01); print comment n power; end; if power >= tpower then do; cell_n=n; omegas=(1-rho**2)*lambda/(lambda+n*j*k*l); if collapse=1 then print 'Power for the BC Interaction'; if collapse=2 then print 'Power for the AC Interaction'; if collapse=3 then print 'Power for the AC Interaction'; comment ='analysis ok'; power=round(power,.01); omegas=round(omegas,.001); print comment omegas cell_n power; n=maxn; end; end; end; quit;