# What is statistical power? And how to conduct power analysis in R?
# R code from https://philippmasur.de/blog/2024-05-28-statistical-power/

# Load package
library(tidyverse)

# Set seed to ensure reproducibility
set.seed(5)

# Simulate true effect in the population
n <- 1e5
control <- rnorm(n, 4, 1)
treatment <- rnorm(n, 5, 1)
pop <- tibble(condition = rep(c("control", "treatment"), each = n),
              outcome = c(control, treatment))
head(pop)

p1 <- ggplot(pop, aes(x = outcome, fill = condition)) +
  geom_density(alpha = .3) +
  theme_void() +
  theme(legend.position = "bottom")
p1

set.seed(11)
exp <- sample_n(pop, size = 10)
head(exp)

p1 +
  geom_point(aes(x = outcome, y = 0, fill = condition),
             data = exp, alpha = .8, shape = 21, size = 5, color = "black")

t.test(outcome ~ condition, exp)

set.seed(11)
exp2 <- sample_n(pop, size = 100)

p1 +
  geom_point(aes(x = outcome, y = 0, fill = condition),
             data = exp2, alpha = .8, shape = 21, size = 5, color = "black")

t.test(outcome ~ condition, exp2)

## n <- 50      # participants per group
## d <- 0.3     # true effect size (Cohen's d)
## 
## p_values <- replicate(1000, {
##   control   <- rnorm(n, mean = 0, sd = 1)
##   treatment <- rnorm(n, mean = d, sd = 1)
##   t.test(treatment, control)$p.value
## })
## 
## mean(p_values < .05)                          # simulated power
## power.t.test(n = n, delta = d, sd = 1)$power  # analytical power

# Load package
library(pwr)

# Computing Cohen's d based on the results of our first experiment
cohens_d <- (4.79-3.99)/1.43

# Estimate power
pwr.t.test(d = cohens_d,                # Effect size
           n = 5,                       # Sample size per group
           sig.level = .05,             # Alpha error
           alternative = "two.sided")   # Type of test

# Computing Cohen's d based on the results of our first experiment
cohens_d <- (5.02-4.05)/1.13

# Post-hoc power analysis
pwr.t.test(d = cohens_d,
           n = 50,
           sig.level = .05,
           alternative = "two.sided")

# Estimate necessary sample size
pwr.r.test(r = .1,                      # Assumed effect size
           sig.level = .05,             # alpha error
           power = .95,                 # Power (1 - beta error)
           alternative = "two.sided")   # Type of test

# Set power range
power <- seq(.50, .95, by = .05)

# Create function to extract sample size from pwr.r.test
extract_n <- function(x) {
  pwr.res <- pwr.r.test(r = .1,                      
                        sig.level = .05,            
                        power = x,                 
                        alternative = "two.sided")
  pwr.res$n
}

# Combine in table
pwr.table <- tibble(power = power,
       sample_size = power |> 
  map_dbl(function(x) extract_n(x)))
pwr.table

# Plot sample curve
pwr.table |> 
  ggplot(aes(x = sample_size, y = power)) +
  geom_line(color = "red") +
  geom_hline(yintercept = .8, linetype = "dashed", color = "grey") +
  theme_classic() +
  labs(x = "Sample Size", y = "Power")

power <- seq(.60, .95, by = .05)
effect <- seq(.1, .5, by = .05)
grid <- expand_grid(power, effect)

# Extract sample sizes
extract_param <- function(x, y) {
  pwr.res <- pwr.r.test(r = y,                      
                        sig.level = .05,            
                        power = x,                 
                        alternative = "two.sided")
  pwr.res$n
}

# Plot curves
tibble(grid, 
       n = pmap_dbl(list(grid$power, grid$effect), extract_param)) |> 
  ggplot(aes(x = n, y = power, group = effect, color = effect)) +
  geom_hline(yintercept = .8, linetype = "dashed", color = "grey") +
  geom_line() +
  theme_classic() +
  labs(x = "Sample Size", y = "Power")

# Load package
library(MonteCarlo)

# Ensure reproducibility
set.seed(42)

# Define parameters
n <- 500
means <- c(2.5, 2.75, 3, 4)  # a1b1, a2b1, a1b2, a2b2
sd <- 1

# Simulate data
dv <- rnorm(n, mean = means , sd = sd) %>% round(0)
iv1 <- rep(c("a1", "a2"),each = 1, n/2)
iv2 <- rep(c("b1", "b2"), each = 2, n/4)
d <-data.frame(dv, iv1, iv2)

# Check data simulation
d %>%
  group_by(iv1, iv2) %>%
  summarise(m = mean(dv)) %>%
  ggplot(aes(x = factor(iv1), y = m, color = iv2, group = iv2)) +
  geom_point(size = 4) +
  geom_line(linewidth = 1) +
  ylim(1, 5) +
  theme_classic() +
  labs(y = "DV", x = "IV 1",
       color = "IV 2")

# Simulation function
sim_func <- function(n = 500, 
                     means = c(2.5, 2.75, 3, 4),
                     sd = 1) {
  
  # Simulate data
  dv <- round(rnorm(n, mean = means , sd = sd),0)
  iv1 <- rep(c("a1", "a2"),each = 1, n/2)
  iv2 <- rep(c("b1", "b2"), each = 2, n/4)
  d <-data.frame(dv, iv1, iv2)
  
  # Fit models
  fit1 <- lm(dv ~ iv1, d)     # Main effet of IV1
  fit2 <- lm(dv ~ iv2, d)     # Main effect of IV2
  fit3 <- lm(dv ~ iv1*iv2, d) # Interaction between both
  
  # Extract p-values and compute significance
  p_1 <- summary(fit1)$coef[2,4]
  sig_1 <- ifelse(p_1 < .05, TRUE, FALSE)
  p_2 <- summary(fit2)$coef[2,4]
  sig_2 <- ifelse(p_2 < .05, TRUE, FALSE)
  p_3 <- summary(fit3)$coef[4,4]
  sig_3 <- ifelse(p_3 < .05, TRUE, FALSE)
  

  # return values as list
  return(list("p_1" = p_1,
              "sig_1" = sig_1,
              "p_2" = p_2,
              "sig_2" = sig_2,
              "p_3" = p_3,
              "sig_3" = sig_3))
}

# Parameters that should vary across simulations
n_grid <- seq(100, 600, 100)
sd_grid <- c(1, 1.5, 2)

# Collect simulation parameters in list
(param_list <- list("n" = n_grid,
                    "sd" = sd_grid))

result <- MonteCarlo(func = sim_func,            # pass test function
                     nrep = 1000,                 # number of tests
                     ncpus = 4,                   # number of cores to be used
                     param_list = param_list, )   # provide parameters

df <- MakeFrame(result)
head(df)

df %>%
  dplyr::select(-contains("p")) %>%
  gather(key, value, -n, -sd) %>%
  group_by(n, sd, key) %>%
  summarize(power = sum(value)/length(value)*100) %>%
  ggplot(aes(x = n, y = power, color = key)) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 80, linetype = "dashed") +
  geom_hline(yintercept = 95, linetype = "dashed") +
  facet_wrap(~sd) +
  theme_bw() +
  labs(title = "Power Curves",
       x = "sample size (n)", 
       y = "power (1-beta)", 
       color = "type of effect",
       caption = "Note: facets represent different standard deviations")