## ------------------------------------------------------------------------
# population SDs
sigA <- 2
sigB <- 4

# sample sizes
n <- 5
m <- 15

# standard deviation of Xbar - Ybar
se <- sqrt(sigA^2/n + sigB^2/m)

# difference between population means
delta <- 3

# critical value for alpha = 0.05
C <- qnorm(0.975)

# the power 
1 - pnorm(C - delta/se) + pnorm(-C - delta/se)

## ------------------------------------------------------------------------
# same parameters as above, but we don't know sigA and sigB

# critical value for alpha = 0.05
C <- qt(0.975, n+m-2)

# the power 
1 - pt(C, n+m-2, delta/se) + pt(-C, n+m-2, delta/se)

## ------------------------------------------------------------------------
# A
power.t.test(n=10, delta=5, sd=10)
# B
power.t.test(delta=5, sd=10, power=0.8)
# C
power.t.test(delta=5, sd=10, power=0.8, alternative="one.sided")

## ------------------------------------------------------------------------
# A function to simulate normal data and get the p-value from the t-test.

# n,m = sample sizes
# sigA,sigB = population SDs
# delta = difference between population means

sim <- function(n, m, sigA, sigB, delta){
  x <- rnorm(n, delta, sigA)
  y <- rnorm(m, 0, sigB)
  t.test(x,y)$p.value
}

# For n=5, m=10, sigA = 1, sigB = 2, and delta = 0, 0.5, 1.0, ..., 2.5, run 1000 simulations and get the p-value:

delta <- c(0.0, 0.5, 1.0, 1.5, 2.0, 2.5)
p <- matrix(nrow=1000,ncol=length(delta))
for(i in 1:length(delta)){
  for(j in 1:nrow(p)){
    p[j,i] <- sim(5, 10, 1, 2, delta[i])
  }
}

# Estimate the power
pow <- apply(p, 2, function(a) mean(a<0.05))
print(pow)
plot(delta,100*pow,xlab=expression(Delta),ylab="power [ % ]",ylim=c(0,100),pch=20,col="blue",cex=1.5)
abline(h=seq(0,100,10),lty="dotted",col="lightgrey")