/* (Programmed by Phil Chapman; revised 9-7-95; send comments or

  corrections to: pchapman@lamar.colostate.edu)

  This program calculates power for the F-test in a randomized complete block

  design with t treatments having true means mu1 to mut, and n blocks.

  If t=2, it is equivalent to a paired (two-sided) t-test.

  Power is the probability of rejecting the null hypothesis for a given set of

  values of the alternative hypothesis.  Power, as we define it,

  (unlike some programs) includes the probability that mean1

  will be concluded to be larger than mean2, when in fact mean1 is smaller than

  mean2.  The probability that such an erroneous conclusion will occur becomes

  negligible as the true difference between the two means becomes large.

  The estimate of experimental error standard deviation within blocks is sigma.

  The size of the test is alpha.

  Power is printed and plotted versus n for values n1 to n2 in increments

  of n3.                        */



data powerrcb;

* user inputs;

t=4;          * t=number of treatments;

              * The following calculates the corrected sum of squares of the true

                treatment means.  It is sometimes called the total effect size.

                There should be t treatment means in the parentheses, separated by commas.;

a2=css(4,5,6,7);

sigma=3;       *sigma= the error standard deviation;

n1=2;          *n1 is the smallest value for n considered;

n2=20;         *n2 is the largest value for n considered;

n3=1;          *n3 is the increment in n;

alpha=0.05;    *alpha is the size (type I error prob) of the test;

* end user inputs;

dfnum=t-1;

do n=n1 to n2 by n3;

dfdenom=(t-1)*(n-1);

fcrit=finv(1-alpha,dfnum,dfdenom);  *fcrit is the F value you would get from the F-table;

ncp=n*a2/sigma**2;       *ncp  is noncentrality parameter (it is divided by 2 in most books;

power=1-probf(fcrit,dfnum,dfdenom,ncp);

output;

end;

proc print;

    var sigma a2 n power;

proc plot;

  plot power*n;

run;