Formulas to calculate power of a test. What am I doing wrong?

Question

In Jay L. Devore 'Probability and Statistics for Engineering and the Sciences', he presents the following formula to retrieve β, when we consider a two-tailed test and the population mean is known.

$\beta_{(\mu')} = \Phi({z_{\alpha/2}+\frac{\mu_0 - \mu'}{\sigma/\sqrt{n}}}) - \Phi({-z_{\alpha/2}+\frac{\mu_0 - \mu'}{\sigma/\sqrt{n}}})$

Therefore I tried to apply this to a mock-up data:

Desired $\alpha = 0.05$

$\mu_0 = 14$

$\mu' = 15$

$n = 4$

$\sigma = 1.496$

From this I get:

$\beta = \Phi(1.96 + \frac{-1}{0.374}) - \Phi(- 1.96 + \frac{-1}{0.374})$

$\beta = \Phi(-0.714) - \Phi(-4.634)$

$β = 0.24$

However, when I try to get power from R language:

> power.t.test(n = 4, delta = 1, sd = 1.496, 
               sig.level = 0.05, power = NULL, 
               type = "one.sample", alternative = "two.sided")

I get power = $0.1578781$ ($\beta = 0.84$).

This is a huge difference and it's making me nuts. What am I doing wrong, and why both formulas seem to be the same thing, though they are obviously not?

Answer 1

There are several mistakes here, first of all, since $n=4$, you should divide $\sigma$ by $\sqrt{4}=2$ instead. So the correct type II rate is

> pnorm(qnorm(0.975) - 2 / 1.496) - pnorm(qnorm(0.025) - 2 / 1.496)
[1] 0.7328904

To check the result in R, you may call pwr.norm.test from package pwr:

> pwr.norm.test(d = 1 / 1.496, n = 4, sig.level = 0.05)

     Mean power calculation for normal distribution with known variance 

              d = 0.6684492
              n = 4
      sig.level = 0.05
          power = 0.2671096
    alternative = two.sided

On the other hand, power.t.test computes the power of the $t$-test, where population $\sigma$ is unknown and must be estimated. In this case, the power function is given by $$F\left(t_{n-1,\alpha_0/2}\right)+1-F\left(-t_{n-1,\alpha_0/2}\right),$$ where $F$ denotes cdf of the non-central $t$-distribution with df $n-1$ and non-centrality $$\psi=\frac{\mu-\mu_0}{\sigma/\sqrt{n}}.$$ In R, the power can be computed in the following ways:

> pt(qt(0.025, 3), 3, ncp = 2 / 1.496) + pt(qt(0.975, 3), 3, ncp = 2 / 1.496, lower.tail = F)
[1] 0.1592488
> power.t.test(n = 4, delta = 1, sd = 1.496, sig.level = 0.05, type = "one.sample" , strict = T)

     One-sample t test power calculation 

              n = 4
          delta = 1
             sd = 1.496
      sig.level = 0.05
          power = 0.1592488
    alternative = two.sided

The $0.1578781$ value you previously got only accounts for one tail (without setting strict = TRUE).