Power analysis

library(ggplot2)
library(tidyverse)

## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr

## Conflicts with tidy packages ----------------------------------------------

## filter(): dplyr, stats
## lag():    dplyr, stats

creating numbers for effect size, sample size, significance level, and sd

effect <- c(15.01, 15.05, 15.1, 15.2, 15.5)
sample <- c(2, 3, 5, 10, 15)
sig <- c(0.001, 0.01, 0.05, 0.10, 0.15)
sd <- c(0.1, 0.25, 0.5, 1, 2)

Example from the book(p392) effect <- c(15.5) sample <- c(3) sig <- c(0.01) sd <- c(0.25)

MyData = expand.grid(effect = effect, sample = sample, sig=sig, sd=sd)

Now let’s compute the power (IN PROGRESS)

MyData$zstat <- (MyData$effect-15)/(MyData$sd/sqrt(MyData$sample))
MyData$critical <- qnorm(MyData$sig/2, mean = 0, sd = 1) # find the critical value for the wanted significance 

level for a two-tailed test

MyData$threshold_u <- -MyData$critical*MyData$sd/sqrt(MyData$sample)+15 # upper bound 
MyData$threshold_l <- MyData$critical*MyData$sd/sqrt(MyData$sample)+15 # lower bound

MyData$power1 <- 1-pnorm(MyData$threshold_u, mean = MyData$effect, sd = MyData$sd/sqrt(MyData$sample))
# probability to be above the upper threshold; using the sampling distribution
MyData$power2 <- pnorm(MyData$threshold_l, mean = MyData$effect, sd = MyData$sd/sqrt(MyData$sample))
# probability to be below the lower threshold
MyData$power <- MyData$power1 + MyData$power2

summary(MyData$power)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00100 0.02272 0.10850 0.25380 0.30200 1.00000

Illustrating the role of the effect size

ggplot (data=subset(MyData, sample == 3 & sig == 0.01 & sd == 0.25)) + 
          geom_line(mapping = aes(x = effect, y = power))

# Illustrating the role of the sample size
ggplot (data=subset(MyData, effect==15.5 & sig==0.01 & sd==0.25)) +
          geom_line(mapping = aes(x = sample, y = power))

# Illustrating the role of the significance level
ggplot (data=subset(MyData, effect==15.5 & sample==3 & sd==0.25)) +
          geom_line(mapping = aes(x = sig, y = power))

# Illustrating the role of the standard deviation of the population
ggplot (data=subset(MyData, effect==15.5 & sample==3 & sig==0.01)) +
          geom_line(mapping = aes(x = sd, y = power))