# Here's a picture showing the sampling distribution of the mean
# under a null hypothesis that mu=100. (The sample size is 25 and
# the standard deviation is 15.) The vertical lines delineate the
# rejection region.

> par(pin=c(6,4))
> plot(seq(90,120,.1),dnorm(seq(90,120,.1),100,3),type='l',axes=FALSE,
+ xlab="",ylab="")
> abline(h=0)
> 1.96*3
[1] 5.88
> lines(c(105.88,105.88),c(0,dnorm(105.88,100,3)))
> lines(c(100-5.88,100-5.88),c(0,dnorm(100-5.88,100,3)))

# Here is the actual sampling distribution of the mean if
# mu is really 110:

> lines(seq(90,120,.1),dnorm(seq(90,120,.1),110,3))

# There's a tiny chance that sampling from that second distribution
# would result in a mean that was significant in the negative
# direction:

> pnorm(100-5.88, 110, 3)
[1] 6.005335e-08

# And there's a very good chance that it would be significant
# in the positive direction:

> 1-pnorm(105.88, 110, 3)
[1] 0.9151756

# So the power, if mu really is 110, is:

> 1-pnorm(105.88, 110, 3) + pnorm(100-5.88, 110, 3)
[1] 0.9151757

# More realistically, we would use a t test for this situation.
# Here is the noncentrality parameter (true mean - null mean)/(standard error):

> 10/3
[1] 3.333333
> delta <- 10/3
> 
# Here is the sampling distribution of the t statistic if
# the null hypothesis is true. Note that the dt() function call
# has only the value and degrees of freedom:

> plot(seq(-3,7,.1),dt(seq(-3,7,.1), 24),type='l');abline(h=0)

# Here is the sampling distribution if the true mean is really 100.
# Note the addition of the noncentrality parameter, delta, to the
# dt() function call:

> lines(seq(-3,7,.1), dt(seq(-3,7,.1), 24, delta),lty=2)

# Here's the power of the test:

> tcrit <- qt(.975,24)
> tcrit
[1] 2.063899
> pt(-tcrit, 24, delta) + (1-pt(tcrit,24,delta))
[1] 0.892017

# Suppose that instead of this one-sample, post-test-only design,
# we had used a two-sample experimental design, dividing our N between
# the two groups, 12 in one and 13 in the other. We stull assume that
# the true mean is 100. The noncentrality parameter is a bit more
# complicated now:

> delta <- (110-100)/15 * sqrt(12*13/25)
> delta
[1] 1.665333

# And the degrees of freedom change:

> tcrit <- qt(.975, 23)

# But, given those changes, power is still calculated the same way:

> pt(-tcrit, 23, delta) + (1-pt(tcrit,23,delta))
[1] 0.3581033

# How much would power improve if we could afford 25 participants
# in each group?

> delta <- (110-100)/15 * sqrt(25*25/50)
> delta
[1] 2.357023
> tcrit <- qt(.975, 48)
> pt(-tcrit, 23, delta) + (1-pt(tcrit,23,delta))
[1] 0.6381928