################################################################################
# R code to simulate one-sample one-sided power

# Strategy:
  # Do this R times:
  #   draw a sample of size N from the distribution specified by the alternative hypothesis
  #     That is, 25 subjects from a normal distribution with mean 102 and sigma 15
  #   Calculate the mean of our sample
  #   Calculate the associated z-statistic
  #   See whether that z-statistic has a p-value < 0.05 under H0: mu=100
  #      If we reject H0, then set reject = 1, else reject = 0.
  # Finally, the proportion of rejects we observe is the approximate power

n     <-  25;             # sample size of 25
mu0   <- 100;             # null hypothesis mean of 100
mu1   <- 102;             # alternative mean of 102
#mu1   <- 105;             # alternative mean of 105
sigma <-  15;             # standard deviation of normal population

alpha <- 0.05;            # significance level

R     <- 10000;           # Repetitions to draw sample and see whether we reject H0
                          #   The proportion of these that reject H0 is the power

reject <- rep(NA, R);     # allocate a vector of length R with missing values (NA)
                          #   to fill with 0 (fail to reject H0) or 1 (reject H0)
for (i in 1:R) {
  sam <- rnorm(n, mean=mu1, sd=sigma);  # sam is a vector with 25 values

  ybar <- mean(sam);  # Calculate the mean of our sample sam

  z <- (ybar - mu0) / (sigma / sqrt(n));  # z-statistic (assumes we know sigma)
                                          # we could also have calculated the t-statistic, here

  pval <- 1-pnorm(z);   # one-sided right-tail p-value
                        #   pnorm gives the area to the left of z
                        #   therefore, the right-tail area is 1-pnorm(z)

  if (pval < 0.05) {
    reject[i] <- 1;  # 1 for correctly rejecting H0
  } else {
    reject[i] <- 0;  # 0 for incorrectly fail to reject H0
  }

}

power <- mean(reject);  # the average reject (proportion of rejects) is the power
power
# 0.1655  for mu1=102
# 0.5082  for mu1=105



################################################################################
# R code to compute exact two-sample two-sided power
library(pwr) # load the power calculation library

pwr.t.test(n = 15,
    d = c(5,10,15,20,25)/7,
    sig.level = 0.05,
    power = NULL,
    type = "two.sample",
    alternative = "two.sided")

## Results
#
#     Two-sample t test power calculation
#
#              n = 15
#              d = 0.7142857, 1.4285714, 2.1428571, 2.8571429, 3.5714286
#      sig.level = 0.05
#          power = 0.4717438, 0.9652339, 0.9998914, 1.0000000, 1.0000000
#    alternative = two.sided
#
# NOTE: n is number in *each* group



################################################################################
# R code to simulate two-sample two-sided power

# Strategy:
  # Do this R times:
  #   draw a sample of size N from the two hypothesized distributions
  #     That is, 15 subjects from a normal distribution with specified means and sigma 7
  #   Calculate the mean of the two samples
  #   Calculate the associated z-statistic
  #   See whether that z-statistic has a p-value < 0.05 under H0: mu_diff=0
  #      If we reject H0, then set reject = 1, else reject = 0.
  # Finally, the proportion of rejects we observe is the approximate power

n     <-  15;             # sample size of 25
mu1   <- 147;             # null hypothesis English mean
mu2   <- c(142, 137, 132, 127, 122);  # Celt means
sigma <-   7;             # standard deviation of normal population

alpha <- 0.05;            # significance level

R     <- 100000;          # Repetitions to draw sample and see whether we reject H0
                          #   The proportion of these that reject H0 is the power

power <- rep(NA,length(mu2));  # allocate a vector to store the calculated power in

for (j in 1:length(mu2)) {  # do for each value of mu2

  reject <- rep(NA, R);     # allocate a vector of length R with missing values (NA)
                            #   to fill with 0 (fail to reject H0) or 1 (reject H0)

  for (i in 1:R) {
    sam1 <- rnorm(n, mean=mu1   , sd=sigma);  # English sample
    sam2 <- rnorm(n, mean=mu2[j], sd=sigma);  # Celt sample

    ybar1 <- mean(sam1);  # Calculate the mean of our sample sam
    ybar2 <- mean(sam2);  # Calculate the mean of our sample sam

    z <- (ybar2 - ybar1) / (sigma * sqrt(1/n+1/n));  # z-statistic (assumes we know sigma)
                                            # we could also have calculated the t-statistic, here

    pval.Left  <- pnorm(z);    # area under left tail
    pval.Right <- 1-pnorm(z);  # area under right tail
    pval <- 2 * min(pval.Left, pval.Right); # p-value is twice the smaller tail area

    if (pval < 0.05) {
      reject[i] <- 1;  # 1 for correctly rejecting H0
    } else {
      reject[i] <- 0;  # 0 for incorrectly fail to reject H0
    }

  }

  power[j] <- mean(reject);  # the average reject (proportion of rejects) is the power
}

power
# [1] 0.49698 0.97517 0.99994 1.00000 1.00000