/* (Programmed by Phil Chapman; revised 5-7-96; send comments or

  corrections to: pchapman@lamar.ColoState.EDU)

  This program calculates power for the F-tests in a two-factor

  experiment in which the first factor (rows) is completely randomized,

  and the second factor (columns) is a split plot or repeated measure.

  It is appropriate in the situation in which subjects are randomly

  assigned to  r treatments, and then repeatedly measured c at time points.

  It is also appropriate for a split plot design where the whole plot

  factor has been completely randomized.

  The user inputs sigma1, the standard deviation of the true values

  for subjects or whole plots; sigma2, the error standard

  deviation among individual measurements, after accounting for

  all model effects; and the r by c matrix of true treatment means, mu.



   The program calculates power for the F-tests for main effects

   and interactions.

   The program requires that SAS/IML be installed.

*/

proc iml;

*Begin user inputs;

n1=2;         *Smallest number of subj (whole plots) per trt considered;

n2=30;        *Largest number of subj (whole plots) per trt considered;

n3=1;         *Increment in number of subj (whole plots) per trt;

alpha=0.05;   *Type I error rate;

sigma1=sqrt(12);     *Between subjects error standard deviation;

sigma2=4;     *Within subjects error standard deviation;



              *The true means are input as the matrix mu.

               In aan r by c factorial mu will have r rows and

               c columns.  Note that a comma is used to

               separate rows. Add rows if necessary.;

mu={1 2 3 4,

   5 6 7 8};

*End user inputs;

create power var{r c row col rcint n1 n2 n3 alpha sigma1 sigma2};

r=nrow(mu);

c=ncol(mu);

ra=j(c,1,1);

ca=j(1,r,1);

rtot=mu*ra;

ctot=ca*mu;

cor=sum(mu)**2/(r*c);

row=ssq(rtot)/c-cor;

col=ssq(ctot)/r-cor;

rcint=ssq(mu)-row-col-cor;

print row col rcint;

append;

quit;

data power2;

  set power;

sig1sq=sigma1*sigma1;

sig2sq=sigma2*sigma2;

dfrow=r-1;

dfcol=c-1;

dfint=(r-1)*(c-1);

do n=n1 to n2 by n3;

    dfe1=r*(n-1);

    critrow=finv(1-alpha,dfrow,dfe1);

    dfe2=r*(c-1)*(n-1);

    critcol=finv(1-alpha,dfcol,dfe2);

    critint=finv(1-alpha,dfint,dfe2);

    powrow=1-probf(critrow,dfrow,dfe1,n*row/(sig2sq+c*sig1sq));

    powcol=1-probf(critcol,dfcol,dfe2,n*col/sig2sq);

    powint=1-probf(critint,dfint,dfe2,n*rcint/sig2sq);

    output;

end;

proc print;

  var sig1sq sig2sq r c n powrow powcol powint;

run;