# This function reproduces the results of the SAS program, giving power calculations for a 
 # completely randomized design.  The arguments are:
 # alpha is the Type I error of the test
 # ndf is the numerator degrees of freedom
 # partlambda is the part of the noncentrality parameter lambda that is multiplied by the sample size n
 # n1, nlast, and ninc are the starting, ending, and increment values for the sample size   
 

powercr <- function(alpha,ndf,partlambda,n1,nlast,ninc)
     {
        nnumber <- round( (nlast -n1)/ninc)
        n       <- numeric(nnumber)
        power   <- numeric(nnumber)
        for (i in seq(n1,nlast,by=ninc))
           {
              n[i]        <-  i
              ddf         <-  (ndf +1)*(i-1)
              fcr         <-  qf(1-alpha,ndf,ddf)
              lambda      <-  i*partlambda
              power[i]    <-  1 - pf(fcr,ndf,ddf,lambda) 
            }
        print(cbind(n,power))
        plot(n,power,type="b")
      }


powercr(.05,3,2,2,50,1)           # up to sample size 50

powercr(.05,3,2,2,20,1)           # up to sample size 20