# This function reproduces the results of the SAS program, giving power calculations for a 
 # completely randomized design with random effects.  The arguments are:
 # alpha is the Type I error of the test
 # ngrp is the number of factor levels (t) in the completely randomized design
 # rat is equal to:  sigma_a^2/sigma_e^2
 # r1, rlast, and rinc are the starting, ending, and increment values for the sample size   
 

powercrRandom <- function(ngrp,rat,alpha,r1,rlast,rinc)
     {
        rnumber <- round( (rlast -r1)/rinc)
        r       <- numeric(rnumber)
        power   <- numeric(rnumber)
        for (i in seq(r1,rlast,by=rinc))
           {
              r[i]        <-  i
              ndf         <-  ngrp -1
              ddf         <-  (ngrp)*(i-1)
              fcr         <-  qf(1-alpha,ndf,ddf)
              lambda2     <-  1 + i*rat
              power[i]    <-  1 - pf(fcr/lambda2,ndf,ddf) 
            }
        print(cbind(r,power))
        plot(r,power,type="b")
      }


powercrRandom(5,.8225,.05,2,15,1)           # up to sample size 15