Monday 18 March 2013

Power calculations using R

Calculating the power of statistical tests

You can calculate the power of a statistical test using R. This is the probability of rejecting the null hypothesis, given that the null hypothesis is false (ie. the situation where you want to reject the null hypothesis).

One-sided Z-test for paired data
We can calculate the power of a Z-test for paired data by assuming the population standard deviation is equal to the sample standard deviation. We can do this in R using the following function:
> calcPowerOfOneSidedZtestForPairedData <- function(n, sd, mu, mu_0, level)
   {
         tstatistic <- (mu - mu_0)/(sd/sqrt(n))
         quantile <- qnorm(1-level)
         power <- pnorm(quantile - tstatistic, lower.tail=FALSE)
         print(paste("power=",power))
   }
For example, if we have 10 paired data points, the sample standard deviation is 1.789, the true mean (mu) is 0.5, and the mean under the null hypothesis (mu_0) is 0, then using a 5% significance level, the power can be calculated as:
> calcPowerOfOneSidedZtestForPairedData(10, 1.789, 0.5, 0, 0.05)
[1] "power= 0.223315962259717"
That is, the power is 0.223, rounded to three decimal places.
What about if we use a 1% significance level?
> calcPowerOfOneSidedZtestForPairedData(10, 1.789, 0.5, 0, 0.01)
[1] "power= 0.0745755634235717"
That is, the probability of rejecting the null hypothesis is lower when we use a 1% significance level than when we use a 5% significance level. That is, the power of the test is smaller.

If you want to calculate the sample size required to get a certain power for a one-sided test, you can use the function calcSampleSizeForOneSidedZtestForPairedData below:
> calcSampleSizeForOneSidedZtestForPairedData <- function(power, sd, mu, mu_0, level)
   {
         quantile <- qnorm(1-level)
         quantile2 <- qnorm(1-power)
         n <- ((sd^2)/((mu - mu_0)^2)) * (quantile - quantile2)^2
         n <- ceiling(n)
         print(paste("n=",n))
   }
For example, to calculate the sample size required to get a power of 0.9 in a one-sided Z-test, if the sample standard deviation is 1.789, and the true mean is 0.5 (mu), and the mean according to the null hypothesis is 0 (mu_0), and we want to use a 5% significance level, we type:
> calcSampleSizeForOneSidedZtestForPairedData(0.9, 1.789, 0.5, 0, 0.05)
 [1] "n= 110"
This tells us that we sample size of at least 110.

One-sided T-test for paired data
If you don't think it's safe to assume that the population standard deviation is equal to the sample standard deviation, then you can calculate the power of a t-test instead. This can be easily done in R using the "power.t.test()" function. For example, if we have 10 paired data points, the sample standard deviation is 1.789, the true mean is 0.5, and the mean according to the null hypothesis is 0, then using a 5% significance level, the power can be calculated as:
> power.t.test(n = 10, delta = 0.5, sd = 1.789, sig.level = 0.05, type=c("one.sample"), alternative=("one.sided"))
   power = 0.20425
This tells us that the power is 0.20425, which is 0.204 rounded to three decimal places.

You can use the same function to calculate the sample size required to get a certain power. For example, to calculate the sample size required to get a power of 0.9 in a one-sided t-test, if the sample standard deviation is 1.789, and the true mean is 0.5, and we want to use a 5% signficance level, we type:
> power.t.test(delta = 0.5, sd = 1.789, power = 0.9,  sig.level = 0.05, type=c("one.sample"), alternative=("one.sided"))
 n = 111.002
This tells us that the sample size must be at least 111.002. Since, the sample size must be whole number, we round this up to 112.

No comments: