# Balíčky
library(WebPower)
library(paramtest)
library(mixedpower)
library(lme4)
library(tidyverse)

# WebPower ----------------------------------------------------
# V zásadě při power analýze operujeme se čtyřmi “veličinami”:
  # hladinou alfa, což je předem stanovená pravděpodobnost chyby I. typu;
  # velikostí účinku;
  # velikostí vzorku;
  # statistickou silou.

# * Korelace --------------------------------------------------

# Můžeme k tomu použít funkci wp.correlation(). Jeden z klíčových 
# argumentů (power, n, r, alpha) musíme nechat prázdný, 
# protože jsme neuvedli power, vypočte se právě ona
wp.correlation(n = 50, r = 0.3, alpha = .05)


wp.correlation(n = c(10, 20, 30), # některé argumenty mohou zahrnovat více hodnot
               r = 0.1,  
               alpha = .05)


wp.correlation(n = seq(50, 300, by = 50),
               r = 0.3, alpha = .05)

# Řada čísel od 0.1 do 0.3 (včetně) po krocích velikosti 0.01
rs <- seq(0.1, 0.3, by = 0.01)

pwr <- rs %>% 
  map(~wp.correlation(power = .8, 
                      r = .x, 
                      alpha = .05)) %>% 
  map_df(~.x[c("r","alpha","power", "n")]) %>% 
  mutate(n = ceiling(n)) %>% 
  print(n = Inf)

pwr %>% 
  ggplot(aes(r, n)) +
  geom_line()

# Co když chceme zjistit power pro různé velikosti vzorku, různé r a 
# různou alfa? Vytvoříme si list relevantních hodnot

power_input <- list(
  n = seq(10, 100, by = 1), # velkost vzorku od 10 do 100
  r = seq(0.1, 0.6, by = 0.1), # velikost korelace od 0.1 po 0.6
  alpha = c(0.05, 0.01) # hladina alfa
)

# Kombinace všech hodnot
power_input <- power_input %>% 
  expand.grid() %>% 
  as_tibble()

# Iterace
power_output <- power_input %>% 
  pmap(wp.correlation)

str(power_output[[1]])

# Extrakce podstatných informací
power_output <- power_output %>% 
  map_df(~.x[c("r", "alpha", "n", "power")])

power_output

# Power plot
power_output %>% 
  ggplot(aes(x = n, y = power, 
             color = as.factor(r))) +
  geom_line(linewidth = 1) +
  facet_wrap(~alpha, labeller = "label_both") +
  geom_hline(yintercept = .8, linetype = 2) +
  labs(color = "r") + 
  scale_x_continuous(breaks = seq(10, 100, by = 10)) +
  scale_y_continuous(breaks = seq(0, 1, by = .1)) +
  theme(legend.position = "top")


# One-Way ANOVA
r2_to_f <- function(r2_full, # Vysvětlený rozptyl po přidání efektu/efektů
                    r2_null = 0, # Původní vysvětlený rozptyl
                    squared = FALSE) { # Chceme čtverec Cohenova f?
  
  # f2 = Přírůstek vysvětleného rozptylu / reziduální rozptyl
  out <- (r2_full - r2_null) / (1 - r2_full)
  
  if (squared == FALSE) {
    out <- sqrt(out)
  }
  
  return(out)
}

# Zpětná konverze
f_to_r2 <- function(f) {
  f^2 / (1 + f^2)
}

cohen_f <- r2_to_f(r2_full = .2)
cohen_f

# Solve for power
wp.anova(k = 4, # Počet skupin
         n = 40, # Celková velikost vzorku
         alpha = .05, # Hladina alfa
         f = cohen_f) # Cohenovo f

# Solve for f and convert to R2
wp.anova(f = NULL, k = 4, n = 100, power = 0.8, alpha = 0.05)
f <- wp.anova(f = NULL, k = 4, n = 100, power = 0.8, alpha = 0.05)$f
f_to_r2(f)

# Solve for n
wp.anova(f = cohen_f, k = 4, n = NULL, power = 0.8, alpha = 0.05)

48/4

# * Lineární regrese ------------------------------
f2 <- r2_to_f(0.25, squared = TRUE)
f2

# Solve for power
wp.regression(n = 50, # velikost vzorku
              p1 = 2, # počet prediktorů
              f2 = f2, # Cohenovo f2
              alpha = 0.01) # Hladina alfa

# Solve for n
f2 <- r2_to_f(0.3, # Celkový vysvětlený rozptyl
              r2_null = 0.25, # Původní vysvětlený rozptyl
              squared = TRUE) # Chceme čtverec Cohenova f
f2

wp.regression(n = NULL, 
              power = .8, # požadovaná statistická síla
              p1 = 3, # Celkový počet prediktorů
              p2 = 2, # Počet prediktorů v původním modelu
              f2 = f2, # Velikost účinku Cohenovo f2
              alpha = 0.01) # Hladina alfa

# Paramtest -------------------------------------------------------
# Obecně musí vytvořená funkce mít následující strukturu

# 1. Kód sloužící ke generování dat.
# 2. Kód sloužící k odhadu modelu (provedení statistického testu) 
#    a uložení tohoto modelu (fit-objektu.
# 3. Kód sloužící k extrakci potřebných informací (statistik) z modelu 
#    (fit-objektu. Obvykle se jedná o klíčové koeficienty modelu, 
#    testovés statistiky, odhady standardních chyb nebo p-hodnoty.
# 4. Stanovení výstupních vektoru. Musí jím být atomický numerický vektoru, 
#    jehož jsou pojmenovány.

t_fun <- function(simNum, # argument sloužící ke stanovení počtu simulací
                  n1, n2, m1, m2, sd1, sd2, # argumenty sloužící pro simulaci vzorků z populace
                  var.equal = FALSE) { # argument používaný funkcí t.test
                                      
  
  x1 <- rnorm(n1, mean = m1, sd = sd1) # simulace dat pro první skupinu
  x2 <- rnorm(n2, mean = m2, sd = sd2) # simulace dat pro druhou skupinu
  
  t <- t.test(x1, x2, var.equal = var.equal)  # provedení t-testu
  
  stat <- unname(t$statistic) # uložení testové statistiky t
  p <- t$p.value # uložení p-hodnoty
  
  return(c(t = stat, p = p)) # output funkce
  # výstupní hodnoty funkce t_fun: named vektor s 
  # testovou statistikou, p-hodnotou a logickou hodnotou informující o tom,
  # zda bychom zamítli nulovou hypotézu
}

pwr_ttest <- run_test(
  t_fun, # použitá funkce
  n.iter = 2000, # počet simulací (kolikrát je funkce spuštěna)
  output = "data.frame", # formát výstupu: list nebo data.frame
  parallel = "snow", # Použití více jader procoru
  ncpus = 7, # Kolik jader použít
  n1 = 20, n2 = 40, m1 = 0, m2 = 0.5, sd1 = 0.5, sd2 = 1.5, # argumenty vložené
  var.equal = FALSE) # do funkce t_fun


res <- results(pwr_ttest) %>% 
  as_tibble()

res

res %>% 
  summarise(power = mean(p <= .05))

pwr_ttest <- grid_search(
  t_fun, # použitá funkce
  params = list(n1 = c(10, 20, 30), 
                n2 = c(40, 50, 60)), # různé velikosti vzorku
  n.iter = 2000, # počet simulací
  output= 'data.frame', # výstupní formát 
  parallel = 'snow', # paralelní zpracování (využití více jader procesoru)
  ncpus = 7, # počet jader procesoru k využití
  m1 = 0, m2 = 0.5, sd1 = 0.5, sd2 = 1.5) # ostatní vstupní argumenty

res <- results(pwr_ttest) %>% 
  as_tibble()

res %>% 
  group_by(n1.test, n2.test) %>% 
  summarise(power = mean(p < 0.05))

res %>% 
  group_by(n1.test, n2.test) %>% 
  summarise(power = mean(p < 0.05)) %>% 
  ggplot(aes(x = factor(n1.test),
             y = power,
             color = factor(n2.test))) +
  geom_point(size = 4) +
  geom_line(aes(group = n2.test)) +
  labs(x = "n1", color = "n2")

# Mixed power ----------------------------------------------
cube <- readRDS("data/cube_data.rds") %>% 
  mutate(word = factor(word) %>% as.integer())

fit <- lmer(y ~ valence + arousal + (1|id) + (1|word),
               data = cube)
summary(fit)

# * Různé velikosti vzorku
pwr1 <- mixedpower(
  model = fit, # použitý mixed model
  data = cube, # Dataset pro fitování modelu
  fixed_effects = c("valence", "arousal"), # Všechny fixní efekty
  simvar = "id", # Náhodný efekt, který chceme variovat
  steps = c(30, 60, 90), # Úrovně zvoleného náhodného efektu
  critical_value = 2, # Kritická hodnota testové statistiky pro zamítnutí H0
  n_sim = 100, # Počet simulací, doporučuje se 1000 či více pro stabilitu výsledků1
  SESOI = c(0, 1, 1), # Velikosti fixních efektů
  databased = FALSE) # Odhad statistické síly na základě efektů vypočtených z dat?

pwr1 %>% 
  pivot_longer(cols = -c(mode, effect), 
               names_to = "sample_size", values_to = "power") %>% 
  ggplot(aes(sample_size, power, group = 1)) +
  geom_point() +
  geom_line() +
  facet_wrap(~effect) +
  expand_limits(y = 0)


# * Různy počet participantů a stimulů 
pwr2 <- R2power(model = fit, 
                data = cube,
                fixed_effects = c("valence", "arousal"),
                simvar = "word", 
                steps = c(10, 20, 30),
                R2var = "id", 
                R2level = 50, # v současné době podpuruje funkce pouze jednu hodnotu druhého náhodného efektu
                critical_value = 2,
                n_sim = 50,
                SESOI = c(0, 1, 1),
                databased = FALSE)
pwr2

pwr2 %>% 
  pivot_longer(cols = -c(mode, effect), 
               names_to = "words", values_to = "power") %>% 
  ggplot(aes(words, power, group = 1)) +
  geom_point() +
  geom_line() +
  facet_wrap(~effect) +
  expand_limits(y = 0) +
  geom_hline(yintercept = .8, linetype = "dashed") +
  scale_y_continuous(breaks = seq(0, 1, by = .2))