Statistical Power

Here we perform a quick experiment to approx statistical power.

library(ggplot2)
set.seed(481)

The data we work with in this quick experiment is generated by the following function.

genData <- function(type = 0, number = 1000) {
    interval.move <- if (type == 0) 
        0 else 0.1
    runif(number, min = -1 + interval.move, max = 1 + interval.move)
}

of the following two forms:

samp <- genData()
summary(samp)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.9990 -0.5030 -0.0137 -0.0068  0.5230  0.9970

ggplot(data = data.frame(samp), aes(x = samp)) + geom_histogram(binwidth = 0.1, 
    position = "identity")

samp <- genData(type = 1)
summary(samp)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -0.895  -0.365   0.102   0.101   0.569   1.100

ggplot(data = data.frame(samp), aes(x = samp)) + geom_histogram(binwidth = 0.1, 
    position = "identity")

Note that the second type of sample comes from a distribution with mean not equal to zero, its mean is \( 0.1 \).

We use a simple test, we reject the hypothesis that the mean is \( 0 \) when the sample mean is greater than \( 0.05 \).

simpleTest <- function(sample) {
    m <- mean(sample)
    if (m > 0.05) 
        FALSE else TRUE
}

runTest <- function(type = 0, number = 1000) {
    sample <- genData(type = type, number = number)
    simpleTest(sample)
}

Now we iteratively apply the simpleTest and consider the failure results.

reps <- 1000
number <- 10  # take 10 samples
test.type <- c(rep(0, 1000), rep(1, 1000))
results <- sapply(X = test.type, FUN = function(x) runTest(x, number))

The following two functions allow to compute the power of these results.

power.table <- function(t) {
    t[1, 2]/(t[1, 2] + t[2, 2])
}

power.tr <- function(testtypes, results) {
    power.table(table(results, testtypes))
}

power.tr(test.type, results)

## [1] 0.575

Now we automate runs of the experiment.

runExper <- function(number = 10, reps = 1000) {
    test.type <- c(rep(0, 1000), rep(1, 1000))
    results <- sapply(X = test.type, FUN = function(x) runTest(x, number))
    list(test.type = test.type, results = results)
}

power.exper <- function(exper) {
    power.tr(exper$test.type, exper$result)
}

uplim <- 500
powers <- sapply(X = seq(from = 1, to = uplim, by = 1), function(n) power.exper(runExper(number = n)))
ggplot(data = data.frame(idx = 1:uplim, powers), aes(x = idx, y = powers)) + 
    geom_point()