# 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
 # ngrp is the number of factor levels (t) in the completely randomized design
 # rat is equal to:  (sum of tau(i)^2)/sigma^2
 # r1, rlast, and rinc are the starting, ending, and increment values for the sample size   
 

powercr <- function(ngrp,rat,alpha,n1,nlast,ninc)
     {
        n       <- seq(n1,nlast,by=ninc)
        power   <- numeric(length(n))
        ndf     <-  ngrp -1
        ddf         <-  ngrp*(n-1)
        fcr         <-  qf(1-alpha,ndf,ddf)
        lambda      <-  n*rat
       # print(c(n,ddf,fcr,lambda))
        power    <-  1 - pf(fcr,ndf,ddf,lambda) 
        print(cbind(n,ddf,lambda,fcr,power))
        plot(n,power,type="b",xlab="sample size per group",ylab="Power")
      }


 # for the text example

powercr(3,8.515,.01,2,10,1)           # up to sample size 15


 # The text discusses using the difference of practical interest, D
 # using this approach, the argument rat above satisfies
 #  rat = D^2/(2*sigma^2)
 

 # ----------------------------------------------------------------------------------------------

 # R has many other tools for power analysis
 # one is the power.anova.test function, included in the base R code
 # For this function you specify the between and within group variances.
 # The between group variances are given by the sample variance of the 
 #  hypothesized true means. 
 # For example, if a one-way CR analysis had true means of 2.82, 3.89, and 3.04
 # and the within group variance is .075
 
var(c(2.82, 3.89, 3.04))
 # = .3193, this is the between group variance

power.anova.test(groups=3,within.var=.075,between.var=.3193,sig.level=.01,n=4)

 # gives the same answers as the function above, or can be easily turned into a power curve:

sampsz  <- seq(2,15,1)
powvals <- sapply(sampsz, function (x) power.anova.test(groups=3,within.var=.075,
                                      between.var=.3193,sig.level=.01,n=x)$power)
plot(sampsz,powvals,type="b")

 # ----------------------------------------------------------------------------------------------

 # there are many R packages for power analysis, one good one is the pwr package