*This program can be used for power calculations for between-subjects main 
effects in completely between-subjects and in split-plot designs.  
Approximations are used to avoid having to specify all variances and covariances
for the repeated measures.  Compound symmetry is assumed. As a result of the approximations, the n will not be exactly 
correct for split plot designs;
options pagesize=65 linesize=65;
data one;
es=.60      ; *Effect size for largest and smallest mean;
alpha=.05   ; *Type I error rate;
j=6         ; *Number of levels for the the between-subjects factor
               for which es is defined. j is limited to 6 or less;
tpower=.80  ; *Target power;
minn= 002   ; *Minimum sample size for calculations;
maxn= 2000  ; *Maximum sample size for calculations;
prodw=1     ; *Product of number of levels of all within-subjects factors.
               Equals p (the number of levels for the the within-subjects 
               factor) when there is just one within-subjects factor.
               Equals 1 when there are no within-subjects factors;
prodb=1      ; *product of number of levels of the other between-subects factors. Equals 1
               for a design with one between-subjects factor design and equals the number 
               of level in the second factor if there are two between-subjects factors;
rho=.0       ; *correlation for within-subjects variables. equals zero if there are 
                no within subjects factors;
prob=1-alpha;
df1=j-1;
if j =3 then do;
  do i = 1 to 2;
    if i = 1 then do;
      config='minimum range and equally spaced                      ';
      mu1 = -es/2;
      mu2 = 0;
      mu3 = es/2;
      delta = (es**2/2)/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x1 cval power mu1 mu2 mu3;
                 do over x1;
                  x1=round(x1,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x2 cval power mu1 mu2 mu3;
                 do over x2;
                  x2=round(x2,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 2 then do;
      config='maximum range, one extreme mean, and two extreme means';
      mu1 = -es/2;
      mu2 = es/2;
      mu3 = es/2;
      delta = ((j**2-1)*es**2/(4*j))/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x3 cval power mu1 mu2 mu3;
                 do over x3;
                  x3=round(x3,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x4 cval power mu1 mu2 mu3;
                 do over x4;
                  x1=round(x4,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
  end;
end;
if j=4 then do;
  do i = 1 to 4;
    if i = 1 then do;
      config='minimum range   ';
      mu1 = -es/2;
      mu2 = 0;
      mu3 = 0;
      mu4 = es/2;
      delta = (es**2/2)/(1+(prodw-1)*rho);
         do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x5 cval power mu1 mu2 mu3 mu4;
                 do over x5;
                  x5=round(x5,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x6 cval power mu1 mu2 mu3 mu4;
                 do over x6;
                  x6=round(x6,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 2 then do;
      config='maximum range and two extreme means  ';
      mu1 = -es/2;
      mu2 = -es/2;
      mu3 =  es/2;
      mu4 =  es/2;
      delta = ((j/4)*es**2)/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x7 cval power mu1 mu2 mu3 mu4;
                 do over x7;
                  x7=round(x7,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x8 cval power mu1 mu2 mu3 mu4;
                 do over x8;
                  x8=round(x8,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 3 then do;
      config='equally spaced  ';
      mu1 = -es/2;
      mu2 = -es/6;
      mu3 =  es/6;
      mu4 =  es/2;
      delta = (j*(j+1)*es**2/(12*(j-1)))/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x9 cval power mu1 mu2 mu3 mu4;
                 do over x9;
                  x9=round(x9,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x10 cval power mu1 mu2 mu3 mu4;
                 do over x10;
                  x10=round(x10,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 4 then do;
      config='one extreme mean';
      mu1=-es/2;
      mu2=-mu1;
      mu3=mu2;
      mu4=mu3;
     delta = ((j-1)*es**2/j)/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x11 cval power mu1 mu2 mu3 mu4;
                 do over x11;
                  x11=round(x11,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x12 cval power mu1 mu2 mu3 mu4;
                 do over x12;
                  x12=round(x12,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
  end;
end;
if j=5 then do;
  do i = 1 to 4;
    if i = 1 then do;
      config='minimum range                             ';
      mu1 = -es/2;
      mu2 = 0;
      mu3 = 0;
      mu4 = 0;
      mu5 = es/2;
      delta = (es**2/2)/(1+(prodw-1)*rho);
         do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x13 cval power mu1 mu2 mu3 mu4 mu5;
                 do over x13;
                  x13=round(x13,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x14 cval power mu1 mu2 mu3 mu4 mu5;
                 do over x14;
                  x14=round(x14,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 2 then do;
      config='maximum range and two extreme means';
      mu1 = -es/2;
      mu2 = -es/2;
      mu3 = -es/2;
      mu4 =  es/2;
      mu5 =  es/2;
      delta = ((j**2-1)*es**2/(4*j))/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x15 cval power mu1 mu2 mu3 mu4 mu5;
                 do over x15;
                  x15=round(x15,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x16 cval power mu1 mu2 mu3 mu4 mu5;
                 do over x16;
                  x16=round(x16,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 3 then do;
      config='equally spaced                            ';
      mu1 = -es/2;
      mu2 = -es/4;
      mu3 = 0;
      mu4 =  es/4;
      mu5 =  es/2;
      delta = (j*(j+1)*es**2/(12*(j-1)))/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x17 cval power mu1 mu2 mu3 mu4 mu5;
                 do over x17;
                  x17=round(x17,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x18 cval power mu1 mu2 mu3 mu4 mu5;
                 do over x18;
                  x18=round(x18,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 4 then do;
      config='one extreme mean                          ';
      mu1=-es/2;
      mu2=-mu1;
      mu3=mu2;
      mu4=mu2;
      mu5=mu2;
      delta = ((j-1)*es**2/j)/(1+(prodw-1)*rho);
           do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x19 cval power mu1 mu2 mu3 mu4 mu5;
                 do over x19;
                  x19=round(x19,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x20 cval power mu1 mu2 mu3 mu4 mu5;
                 do over x20;
                  x20=round(x20,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
  end;
end;
if j=6 then do;
  do i = 1 to 5;
    if i = 1 then do;
      config='minimum range    ';
      mu1 = -es/2;
      mu2 = 0;
      mu3 = 0;
      mu4 = 0;
      mu5 = 0;
      mu6 = es/2;
      delta = (es**2/2)/(1+(prodw-1)*rho);
         do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x21 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x21;
                  x21=round(x21,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x22 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x22;
                  x22=round(x22,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 2 then do;
      config='maximum range    ';
      mu1 = -es/2;
      mu2 = -es/2;
      mu3 = -es/2;
      mu4 =  es/2;
      mu5 =  es/2;
      mu6 =  es/2;
      delta = (j*es**2/4)/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x23 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x23;
                  x23=round(x23,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x24 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x24;
                  x24=round(x24,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 3 then do;
      config='equally spaced   ';
      mu1 = -es/2;
      mu2 = 3*mu1/5;
      mu3 = mu1/5;
      mu6 =  es/2;
      mu4 = mu6/5;
      mu5 = 3*mu6/5;
      delta = (j*(j+1)*es**2/(12*(j-1)))/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x25 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x25;
                  x25=round(x25,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x26 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x26;
                  x26=round(x26,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 4 then do;
      config='one extreme mean ';
      mu1=-es/2;
      mu2=-mu1;
      mu3=mu2;
      mu4=mu2;
      mu5=mu2;
      mu6=mu2;
      delta = ((j-1)*es**2/j)/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x27 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x27;
                  x27=round(x27,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x28 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x28;
                  x28=round(x28,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
    if i = 5 then do;
      config='two extreme means';
      mu1=-es/2;
      mu2= mu1;
      mu3=-mu2;
      mu4=mu3;
      mu5=mu3;
      mu6=mu3;
      delta = (2*(j-2)*es**2/j)/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x29 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x29;
                  x29=round(x29,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x30 cval power mu1 mu2 mu3 mu4 mu5 mu6;
                 do over x30;
                  x30=round(x30,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
     end;
  end;
end;
if j =2 then do;
      mu1 = -es/2;
      mu2 = es/2;
      delta = (es**2/2)/(1+(prodw-1)*rho);
        do n=minn to maxn;
           lambda=n*prodw*prodb*delta;
           df2=j*prodb*(n-1);
           cval=finv(prob,df1,df2);
           power=1-probf(cval,df1,df2,lambda);
             if power lt tpower and n = maxn then do;
               comment = 'raise maxn  ';
               array x31 cval power mu1 mu2 mu3;
                 do over x31;
                  x31=round(x1,.01);
                 end;
               output;
             end;
             if power ge tpower then do;
               comment ='analysis ok';
               array x32 cval power mu1 mu2;
                 do over x32;
                  x32=round(x32,.01);
                 end;
               output;
               n=maxn;
             end;
        end;
    end;
data two;
set one; if j=2;
cell_n=n;
proc sort; by n;
proc print;
Title 'Note that mus are in standard deviation units';
var comment prodw prodb rho es J alpha tpower power mu1 mu2 cell_n lambda cval;

data three;
set one; if j=3;
cell_n=n;
proc sort; by n;
proc print;
Title 'Note that mus are in standard deviation units';
var comment config prodw prodb rho es J alpha tpower power mu1 mu2 mu3 cell_n lambda cval;
data four;
set one; if j=4;
cell_n=n;
proc sort; by n;
proc print;
var comment config prodw prodb rho es J alpha tpower power mu1 mu2 mu3 mu4  cell_n lambda cval;
data five;
set one; if j=5;
cell_n=n;
proc sort; by n;
proc print;
var comment config prodw prodb rho es J alpha tpower power mu1 mu2 mu3 mu4 mu5 cell_n lambda cval;
data six;
set one; if j=6;
cell_n=n;
proc sort; by n;
proc print;
var comment config prodw prodb rho es J alpha tpower power mu1 mu2 mu3 mu4 mu5 mu6 cell_n lambda cval;
proc datasets;
delete one;
quit;
run;