*Power analysis for one-way random effects model with 5 levels;
data power;
alpha=.05;
a=5;
ndf=a-1;
do sig_ratio=0.2 to 3.0 by 0.2; *Ratio of sigma^2_a to sigma^2;
  do n=2 to 25;
    ddf=a*(n-1);
    fcrit=finv(1-alpha,ndf,ddf);
	scaled_fcrit=(1/(n*sig_ratio+1))*fcrit;
    power=1-probf(scaled_fcrit,ndf,ddf);
	output;
  end;
end;
run;
proc print data=power;
run;

*Use this statement if not generating the contour plots in Minitab;
legend1 position=(right middle)
        label=(position=top)
        across=1;
axis1 order=(0 to 25 by 5) minor=(number=4);
axis2 order=(0.0 to 3.0 by 0.5) minor=(number=4);
proc gcontour data=power;
label n="Sample size" sig_ratio="Variance ratio" power="Power";
plot n*sig_ratio=power/levels=(.75 .8 .85 .9 .95)  vaxis=axis1 haxis=axis2 legend=legend1;
run;