################################################################################ #### Demo of Power in Hypothesis Tests #### ################################################################################ #### By Jimin Ding, 02/23/2017 PowerDemo1=function(mu0,mu1,sigma,n,alpha=0.05){ se=sigma/sqrt(n) za=-qnorm(alpha) x=(mu1-4*se):(mu0+4*se) y1=dnorm(x,mean=mu0,sd=se) y2=dnorm(x,mean=mu1,sd=se) y0=rep(0,length(y1)) plot(x,y1,type='l',xlab="",ylab="",xaxt="n",col='red',lwd=2) polygon(c(x,rev(x)),c(y1,y0),col='red',density=20) abline(v=mu0,col='red',lty=2,lwd=2) axis(1,at=mu0,labels=expression(mu[0]),tick=FALSE,line=-0.5) axis(1,at=mu0,labels=mu0,tick=FALSE,line=0.5) readline("This is the null distribution.") xcut=mu0-za*se abline(v=xcut,col='black',lty=1,lwd=3) k=which(xxcut) }else if (alternative=="two.sided"){ xcut=c(mu0-za*se,mu0+za*se) k=which(xxcut[2]) } abline(v=xcut,col='black',lty=1,lwd=3) polygon(c(x[k],rev(x[k])),c(y1[k],y0[k]),col=rgb(1,0,0,alpha=0.7)) if (alternative=="two.sided"){ polygon(c(x[k2],rev(x[k2])),c(y1[k2],y0[k2]),col=rgb(1,0,0,alpha=0.7)) } if (alternative=="less"){ arrows(xcut,-0.00065,xcut,-0.00095,lwd=2,length=0,xpd=TRUE) arrows(xcut,-0.0008,x[1],-0.0008,lwd=2,length=0.1,xpd=TRUE) mtext("Reject Null",side=1,at=xcut-2*se,line=3) arrows(xcut,-0.0008,x[length(x)],-0.0008,lwd=2,length=0.1,xpd=TRUE) mtext("Fail to Reject Null",side=1,at=xcut+2*se,line=3) }else if (alternative=="greater"){ arrows(xcut,-0.00065,xcut,-0.00095,lwd=2,length=0,xpd=TRUE) arrows(xcut,-0.0008,x[1],-0.0008,lwd=2,length=0.1,xpd=TRUE) mtext("Fail to Reject Null",side=1,at=xcut-2*se,line=3) arrows(xcut,-0.0008,x[length(x)],-0.0008,lwd=2,length=0.1,xpd=TRUE) mtext("Reject Null",side=1,at=xcut+2*se,line=3) }else if (alternative=="two.sided"){ arrows(xcut[1],-0.00065,xcut[1],-0.00095,lwd=2,length=0,xpd=TRUE) arrows(xcut[2],-0.00065,xcut[2],-0.00095,lwd=2,length=0,xpd=TRUE) arrows(xcut[1],-0.0008,x[1],-0.0008,lwd=2,length=0.1,xpd=TRUE) mtext("Reject Null",side=1,at=xcut[1]-2*se,line=3) arrows(xcut[2],-0.0008,x[length(x)],-0.0008,lwd=2,length=0.1,xpd=TRUE) mtext("Reject Null",side=1,at=xcut[2]+2*se,line=3) mtext("Fail to Reject Null",side=1,at=mu0,line=3) } abline(v=mu1,col='blue',lty=2,lwd=2) axis(1,at=mu1,labels=mu1,tick=FALSE,line=0.5) lines(x,y2,col='blue',lwd=2) polygon(c(x,rev(x)),c(y2,y0),col='blue',density=20) axis(1,at=mu1,labels=expression(mu[1]),tick=FALSE,line=-0.5) if (alternative!="two.sided"){ polygon(c(x[-k],rev(x[-k])),c(y2[-k],y0[-k]),col=rgb(0,0,1,alpha=0.5)) } else{ polygon(c(x[-c(k,k2)],rev(x[-c(k,k2)])),c(y2[-c(k,k2)],y0[-c(k,k2)]),col=rgb(0,0,1,alpha=0.5)) } polygon(c(x[k],rev(x[k])),c(y2[k],y0[k]),col=rgb(0,1,0,alpha=0.1),border='green') if (alternative=="two.sided"){ polygon(c(x[k2],rev(x[k2])),c(y2[k2],y0[k2]),col=rgb(0,1,0,alpha=0.1),border='green') } legend('topright',c(expression(alpha),expression(beta),"Power"),col=c('red','blue','green'),lty=1,lwd=10) }