power analysis for within-subject, high-repetition experiment

data simulation

library(tidyverse)
library(magrittr)
library(afex)
library(Superpower)

set.seed(555)
SUBJECT_COUNT = 3 # per between-subject condition

#### ground-truth parameters for 3w x 2w x 2b (100 repetitions) 

# get every combination of subjectID, within-subject conditions, and trial repetitions
data = expand.grid(
  subjectID = rnorm(SUBJECT_COUNT), # individual differences
  color = c(-1, 0, 1),
  number = c(-1, 1),
  shape = 1,
  trial_index = 1:100
)

# repeat for the second between subject condition (shape=2)
data = data %>% bind_rows(
expand.grid(
  subjectID = rnorm(SUBJECT_COUNT), # individual differences
  color = c(-1, 0, 1),
  number = c(-1, 1),
  shape = 2,
  trial_index = 1:100
))

#### simulate experiment
data = data %>% 
  # generate responses with color*number interaction with individual differences
  mutate(response = color * number * (15+subjectID*5)) %>% 
  # add trial noise
  mutate(response = response + rnorm(n(), subjectID, 3)) %>% 
  # condition labels
  mutate(
    subjectID=factor(subjectID, labels = paste("Subject", 1:length(unique(subjectID)))), 
    color=factor(color, labels = c("red", "purple", "blue")), 
    number=factor(number, labels = c("one", "three")), 
    shape=factor(shape, labels = c("circle", "square")))

#### preview the data
ggplot(data) +
  aes(x=response, fill=color) +
  geom_density() +
  scale_fill_identity() +
  facet_grid(shape * subjectID ~ number)

#### anova to get partial eta squared (pes)
f_table = afex::aov_ez(
  id="subjectID", 
  dv="response", 
  data, 
  within=c("color", "number"), 
  between="shape",
  anova_table = list(correction="none", es = "pes")) %>% 
  # use pes to calculate cohen's f
  `$`('anova_table') %>% 
  mutate(cohens_f = sqrt(pes/(1-pes))) %>% 
  # why are the row names so hard to access?
  as.data.frame() %>% 
  rownames_to_column("effect") %T>%
  print()

#### power analysis
f_table %>% 
  filter(effect == "color:number") %$%  # fancy pipe
  Superpower::power.ftest(`num Df`, `den Df`, cohens_f)

output.md

 Power Calculation for F-test 

     num_df = 2
     den_df = 8
    cohen_f = 2.888124
alpha_level = 0.05
 beta_level = 1.505091e-05
      power = 99.99849