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

   revisions to: pchapman.colostate.edu)

   This program calculates power for the F-test for a contrast in a one-way

   completely randomized design having true means mu1 to mut, and n observations

   per treatment.  Power, as we define it, includes the probability that a

   contrast will be significantly greater than zero, when it is in fact less

   than zero.   The probability that such an erroneous conclusion will occur

   becomes negligible when the true value of the contrast is large.

  The estimate of experimental error standard

  deviation 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.

  The true means and the contrast coefficients are input as vectors mu and c.

*/

proc iml;

create power var{t n1 n2 n3 sigma alpha a2};

* Begin user inputs;

t=8;          * t=number of treatments;

              * The following is a vector of the true treatment means.

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

mu={1,2,3,4,5,6,7,8};

              * The following is a vector of contrast coefficients.

                There should be t contrast coefficients, separated by commas.;

c={1,1,1,1,-1,-1,-1,-1};

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;

a2=(t(c)*mu)**2/ssq(c);

append;

quit;

data power2;

   set power;

dfnum=1;

do n=n1 to n2 by n3;

dfdenom=t*(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;