This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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