##### Estimating Power Curves




### Preamble
source("http://rfs.kvasaheim.com/math322.R")
library(VGAM)




### Power of t-test
#   Requires data from Normal

## Power when mu wrong

pval = numeric()

for(i in 1:1e4) {
  x = rnorm(15, m=1, s=1)
  pval[i] = t.test(x, m=0)$p.value
}

alpha=0.05
mean(pval<alpha) # This is the power, 0.9477 (for my run)

# Thus, the probability of rejecting the null hyp when mu=1,
# (instead of mu=0) is about 95%.





###
# What about estimating a power curve for
# varying values of mu?

pwr = numeric()
muValues = seq(0,2,length=40) 

for( j in 1:length(muValues) ) {

  pval = numeric()
  thisMU = muValues[j]

  # Est power for this mu
  for(i in 1:1e4) {
    x = rnorm(10, m=thisMU, s=1)
    pval[i] = t.test(x)$p.value  # for hyp of mu = 0
  }
  pwr[j] = mean(pval<alpha) # calculates the power
}

## nb: This took 30 seconds on my desktop






# A plot of power against the mean

par(las=1, family="serif")
par(xaxs="i",yaxs="i")
par(font.lab=2, cex.lab=1.1)
par(mar=c(4,4,0.5,1))

plot.new()
plot.window( xlim=c(0,2), ylim=c(0,1) )

abline(h=0.05, lty=3)
mtext(side=4,at=0.05, text=expression(alpha), line=0.2)
mtext(side=2,at=0.05, text=0.05, line=0.6)
abline(h=0.80, lty=2)

axis(1)
axis(2)
title(xlab=expression(paste("True Mean, ", mu)), line=2.5)
title(ylab="Power", line=2.75)

par(xpd=TRUE)
points(muValues,pwr, pch=21, bg="pink")
par(xpd=FALSE)