Power/Sample Size in Balanced One-Way Random-Effect Model

> # The function power.1rem.test, borrowed from base R function

> # power.anova.test with slight modification, calculates

> # power/sample size for balanced 1-way random effect model.

> # Fetch it (poworem.R) from

> # http://www.stat.wmich.edu/wang/R/codes

> # The following R function can help you get R codes from

> # this location:

> #     "fetch" <-

> #     function(what, where="http://www.stat.wmich.edu/wang/R/codes",

> #              extension="R")

> #     {what <- paste(what, extension, sep=".")

> #      what <- paste(where, what, sep="/")

> #      source(what)

> #     }

> # After this function is created, you can now do

> #     fetch("poworem")

> # Be sure to save the work space following this.

> #

> # Use help(power.anova.test) to see the help documentation

> # for the fixed-effects case (similar to the random-effects

> # case here).

> #

> # Now, you are ready to use the function.

> #

> # With the classnotes example, k (i.e., groups) = 3,

> # σu2 (i.e., between.var) = 50,

> # σ2 (i.e., within.var) = 10, power = .95 and

> # α (i.e., sig.level) = 0.05 (the default), the sample

> # size n can be computed as

> power.1rem.test(groups=3,

+     between.var=50, within.var=10,

+     power=.95)

     Balanced one-way random effect model power calculation 

         groups = 3
              n = 12.53416
    between.var = 50
     within.var = 10
      sig.level = 0.05
          power = 0.95

 NOTE: n is number in each group 

> # Hence the minimum sample size is n = 13 with actual power

> # computed by

> power.1rem.test(groups=3, n=13,

+     between.var=50, within.var=10)

     Balanced one-way random effect model power calculation 

         groups = 3
              n = 13
    between.var = 50
     within.var = 10
      sig.level = 0.05
          power = 0.9518785

 NOTE: n is number in each group

> #

> # The process of finding the desired sample size

> # can be visualized below.

> #

> n <- seq(8,18,.1)

> p <- numeric(length = (len <- length(n)))

> for(i in 1:len) p[i] <- power.1rem.test(3, n[i], 50, 10)$power

> windows(width=4,height=3,pointsize=10)

> plot(p~n, type="l", xlab="sample size", ylab="power",

+      main="Sample Size in Blanced\nOne-Way Random Effect Model",

+      lwd=2)

> abline(h=0.95, lty="dotted", col="red")

> abline(v=13, col="green3")

> axis(1, at=13)

> axis(2, at=.95, label=".95", las=1)

> text(14, .945, adj=0, label=bquote(sigma[u]^2==50))

> text(14, .937, adj=0, label=bquote(sigma^2==10))

> text(14, .929, adj=0, label=bquote(alpha==0.05))

> text(14, .921, adj=0, label="k = 3")

> dev.off()

null device 
          1

(power01.png here)

> #

> # Now find out for the fixed between and within treatment variances

> # as stated (50 and 10, respectively), the minimum number of total

> # sample size N = k×n that will achieve a power of at least

> # 0.95. We begin this by finding sample sizes for k = 2, 3, 4, 5,

> # 6, and 7.

> #

> n <- numeric(6)

> for(i in 1:6)

+   n[i] <- power.1rem.test(i+1, between=50, within=10, power=.95)$n

> names(n) <- as.character(2:7)

> # The command below gives the required sample sizes.

> ceiling(n)

  2   3   4   5   6   7 
197  13   6   4   3   3 

> # The command below gives the required total sample sizes.

> (2:7) * ceiling(n)

  2   3   4   5   6   7 
394  39  24  20  18  21

> # It appears that for k = 6 groups, and n = 3 replications,

> # the minimum total sample size is attained with a value of

> # N = 18.

> #

> # Note that the actual powers at the described numbers of

> # groups and required replications are computed below:

> #

> power.1rem.test(2:7,ceiling(n),50,10)$power

        2         3         4         5         6         7 
0.9501080 0.9518785 0.9591045 0.9621973 0.9589242 0.9783840

> #

> # Now plot the power curve for k = 6, n = 3, σ2 = 10

> # for various values of σu2:

> su <- 25:100

> p <- power.1rem.test(6, 3, su, 10)$power

> windows(width=4, height=3, pointsize=10)

> plot(p~su, type="l", lwd=2,

+   xlab="", ylab="power",

+   main="Power Curve for Balanced\nOne-Way Random Effect Model")

> title(xlab=bquote(sigma[u]^2))

> abline(v=50, lty="dotted", col="green3")

> # The command below shows the power at σu2=50.

> p[su==50]

[1] 0.9589242

> abline(h=p[su==50], lty="dotted", col="red")

> axis(1, at=50)

> text(85, p[su==50], pos=3, offset=0.2, col="red",

+      label=bquote("power"==.(round(p[su==50],4))))

> text(60, .93, adj=0,

+   label=bquote(sigma^2==10))

> text(60, .91, adj=0, label="k = 6")

> text(60, .89, adj=0, label="n = 3")

> dev.off()

null device 
          1

(power02.png here)

> #

> # Now let's graph a contour map of the power function at

> # various values of σ2 and σu2

> #

> # First define grids by giving various values of these two variance

> # components.

> su <- seq(40,60,length=41)

> s <- seq(5, 15, length=41)

> #

> # The command below produces a power table at the grid points.

> p <- outer(s, su,

+   function(x,y) power.1rem.test(6,3,y,x)$power)

> #

> # Now make the graph.

> windows(width=4,height=3,pointsize=10)

> par(pty="s")

> contour(s, su, p, main="Contour Map of Power",

+   xlab=bquote(sigma^2), ylab=bquote(sigma[u]^2))

> abline(v=10,h=50,lty="dotted",col="red")

> par(xpd=NA)

> text(par('usr')[2],60,label="k = 6\nn = 3",

+      adj=-0.2)

> dev.off()

null device 
          1

(power03.png here)