Skip to content

Instantly share code, notes, and snippets.

@dgrapov
Last active December 11, 2015 05:48
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save dgrapov/4554191 to your computer and use it in GitHub Desktop.
#purpose: model the relationship between effects size, sample size and power
# holding the significance level constant (p = 0.05) for three of the most common statistical tests:
# 1) t-Test (two-sample)
# then visualize the results for medium small, medium and large effect sizes (0.2,0.5,0.8 ; http://www.uccs.edu/lbecker/effect-size.html)
#Need pwr package
if(!require(pwr)){install.packages("pwr");library("pwr")}
# t-TEST
#---------------------------------
d<-seq(.1,2,by=.1) # effect sizes
n<-1:150 # sample sizes
t.test.power.effect<-as.data.frame(do.call("cbind",lapply(1:length(d),function(i)
{
sapply(1:length(n),function(j)
{
power.t.test(n=n[j],d=d[i],sig.level=0.05,power=NULL,type= "two.sample")$power
})
})))
t.test.power.effect[is.na(t.test.power.effect)]<-0 # some powesr couldn't be calculated, set these to zero
colnames(t.test.power.effect)<-paste (d,"effect size")
#plot results using base
#------------------------------------------------
obj<-t.test.power.effect # object to plot
imp<-c(2,5,8) # cuts
cuts<-c("small","medium","large") # based on cohen 1988
cols<-1:ncol(obj)
color<-rainbow(length(cols), alpha=.5) # colors
lwd=5 # line thickness
lty<-rep(1,length(color))
lty[imp]<-c(2:(length(imp)+1))
#highligh important locations
imp<-c(2,5,8) # cuts
cuts<-c("small","medium","large") # based on cohen 1988
color[imp]<-c("black")
wording<-d
wording[imp]<-cuts
par(fig=c(0,.8,0,1),new=TRUE)
#initialize plot
plot(1,type="n",frame.plot=FALSE,xlab="sample size",ylab="power",xlim=c(1,150),ylim=c(0,1),main="t-Test", axes = FALSE)
#add custom axis and grid
abline(v=seq(0,150,by=10),col = "lightgray", lty = "dotted")
abline(h=seq(0,1,by=.05),col = "lightgray", lty = "dotted")
axis(1,seq(0,150,by=10))
axis(2,seq(0,1,by=.05))
#plot lines
for(i in 1:length(cols)){lines(1:150,obj[,cols[i]],col=color[i],lwd=lwd,lty=lty[i])}
#legend
par(fig=c(.65,1,0,1),new=TRUE)
plot.new()
legend("top",legend=wording,col=color,lwd=3,lty=lty,title="Effect Size",bty="n")
#plot using ggplot2
#------------------------------------------------
#plot results using ggplot2
library(ggplot2);library(reshape)
x11()
obj<-cbind(size=1:150,t.test.power.effect) #flip object for melting
melted<-cbind(melt(obj, id="size"),effect=rep(d,each=150)) # melt and bind with effctfor mapping
ggplot(data=melted, aes(x=melted$size, y=melted$value, color=as.factor(melted$effect))) + geom_line(size=2,alpha=.5) +
ylab("power") + xlab("sample size") + ggtitle("t-Test")+theme_minimal()
# wow ggplot2 is amazing in its brevity
# need to tweak legend and lty, but otherwise very similiar
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment