Sample size determination in R, #2

Balanced One-Way Fixed Effects Model

Given maximal treatment means difference (effect size) Δ=9, experimental error σ=3, significance level α=0.01, and targeted power γ=0.90, find the common group sample size n.

> source("http://www.stat.wmich.edu/wang/664/egs/power.R")  # do this only once

> pow(6,(delta<−9)/(sigma<−3))   # power when n=6

[1] 0.8976995

> (pow(7,delta/sigma)−>actual)     # power when n=7

[1] 0.9574194

> sapply(as.list(5:8),function(n)pow(n,delta/sigma))  # powers when n=5:8

[1] 0.7784404 0.8976995 0.9574194 0.9837236

# Next command finds required sample size for power=0.90

> uniroot(function(x) sapply(x,function(n)pow(n,delta/sigma))-0.9,

+         lower=5,upper=8)$root

[1] 6.027413

>

# Next command finds powers at Delta=seq(0,10,.1) when n = 7

> y <− sapply(as.list(x<−seq(0,10,.1)),function(d)pow(7,d/sigma))

# Now, produce power curve

> windows(width=5, height=5, pointsize=11)

> plot(y~x,type="l",xlab=bquote(Delta),ylab="Power", lwd=2,

+      main="Power Curve\nn=7")

> abline(h=.9,col=rgb(.75,.75,.75,.7))

> abline(v=9,col=rgb(0,1,1,.7))

> abline(h=actual, col=rgb(.93,.56,.56,0.7))

> text(2,.9,pos=1,lab="targeted power = 0.9",col=rgb(.75,.75,.75,.5))

> text(2,actual,pos=3, col=rgb(.93,.56,.56,0.5),

+      lab=paste("actual power",round(actual,3),sep=" = "))

> text(9, par('usr')[3],pos=1,label="9",xpd=NA)

> dev.off()

null device 
          1 
(Rpower01.png here)

>

> # Now, try to produce several power curves at various group sizes

> f <− function(x, gsize, sigma=3, alpha=0.01)

+        sapply(as.list(x), function(d)pow(gsize,d/sigma,alpha))

> line.types <− line.cols <− 1:4

> gs <− c(7, 5,6,8);  line.widths <− c(2, rep.int(1,3))

> # submit the following commands one at a time to see how a graph

> # be constructed, element by element

> plot.new()   # start a new graph

> plot.window(xlim=c(0,10), ylim=0:1)   # set up coordinates

> abline(v=0:10, col=rgb(.9,.9,.9,0.8))  # vertical grids

> abline(h=0.01, col=rgb(.7,.7,.7,.8))   # significance level

> text(8,0.01,pos=3,col=rgb(.7,.7,.7,.8),

+      label=expression(paste("significance level ",alpha==0.01)))

> abline(h=.9,col=rgb(.75,.75,.75,.7))   # targeted power

> text(2,0.9,pos=1,col=rgb(.75,.75,.75,.7),

+      label=expression(paste("targeted power ",gamma==0.90)))

> abline(h=actual,

+        col=rgb(.93,.56,.56,0.7))  # actual power at n=7, delta=9

> text(2,actual,pos=3,col=rgb(.93,.56,.56,0.7),

+      label=bquote(paste("actual power ",gamma*"*"==.(round(actual,3)))))

> abline(v=9, col=rgb(.75,.75,.75,.7))  # effect size, Delta=9

> text(9, par('usr')[3], pos=1, col=rgb(.75,.75,.75,.7), xpd=NA,

+      label=bquote(Delta==9))

> axis(1);  axis(2, las=1)   # two axes

> title(xlab="Effect Size", ylab="Power")   # axes' titles

> title(main="Power Curves at Various Group Sample Sizes")  # main title

> title(sub="last paragraph, page C.13", adj=0)   # subtitle

> mtext(side=3,line=0, text=expression(paste(Delta==9,", ", sigma==3,

+       ", ", gamma==0.90)) )  # margin text at side 3 (top)

> box()   # box plot region

> # now, add power curves one at a time

> for (i in 1:4)

+   plot(function(x)f(x,gs[i]), from=0, to=10, add=T,

+          lty=line.types[i], col=line.cols[i], lwd=line.widths[i])

> legend(0,0.6, legend=as.character(gs), title="Group Sample Size", ncol=2,

+               lty=line.types, lwd=line.widths, col=line.cols,

+               bg=rgb(.9,.93,.93,.8))  # add a legend

> dev.off()

null device 
          1 
(Rpower02.png here)

>

> rm(x,y,actual,i)   # cleanup