# Balíčky a custom funkce -------------------------------------------------
library(WebPower)
library(paramtest)
library(tidyverse)


pwr_analysis <- function(input, fn, ...) {
  
  if (!is.data.frame(input)) {
    input <- input %>% 
      purrr::cross_df()
  }
  
  out_list <- purrr::pmap(input, fn, ...)
  
  out_df <- out_list %>% 
    purrr::map_df(~.x[names(.x)])
  
  out_df
}

pwr_plot  <- function(data, x = "n", y = "power", 
                      color = NULL, facet_wrap = NULL, facet_grid = NULL,
                      line_width = 1,
                      hline = NULL, vline = NULL) {
  
  # Převést sloupce color na factor, pokud ještě není factor
  if (!is.null(color) && !is.factor(data[[color]])) {
    data[[color]] <- factor(data[[color]])
  }
  
  # Vytvořit základní kostru grafu
  p <- data %>% 
    ggplot(aes_string(x = x, y = y)) +
    theme(legend.position = "top") +
    expand_limits(y = 0)
  
  # Přidat křivku. Pokud je color vyplněno, přidat křivky různé barvy
  if (is.null(color)) {
    p <- p + geom_line(aes(group = 1), size = line_width)
  } else {
    p <- p + 
      geom_line(aes_string(color = color, group = color), size = line_width)
  }
  
  # Přidat horizontální nebo vertikální čáru
  if (!is.null(hline)) {
    p <- p + geom_hline(yintercept = hline,
                        linetype = 2)
  }
  
  if (!is.null(vline)) {
    p <- p + geom_hline(yintercept = vline,
                        linetype = 2)
  }
  
  # Přidat fazety, pokud jsou vyplněny
  if (!is.null(facet_wrap)) {
    p <- p + facet_wrap(vars(.data[[facet_wrap]]), 
                        labeller = "label_both")
  } else if (!is.null(facet_grid)) {
    p <- p + facet_grid(rows = vars(.data[[facet_grid[1]]]),
                        cols = vars(.data[[facet_grid[2]]]),
                        labeller = "label_both")
  }
  
  return(p)
}

# Korelace ---------------------------------------------------

# Spočítejte statistickou sílu pro oboustranné testy tří korelačních koeficientů 
# r = 0,1, 0,2 a 0,3 při hladině alfa = 0,05 a velikosti vzorku od 30 do 1000 
# a vytvořte vhodné grafy ukazující, jak statistická síla těchto 
# tří korelačních koeficientů závisí na velikosti vzorku. Klidně k tomu 
# použijte funkce pwr_analysis() a pwr_plot().
pwr_cor <- pwr_analysis(
  list(r = c(0.1, 0.2, 0.3),
       n = 30:1000),
  wp.correlation
)

pwr_plot(pwr_cor,
         x = "n",
         y = "power",
         color = "r") +
  scale_x_continuous(breaks = seq(0, 1000, by = 50))

pwr_cor %>% 
  group_by(r) %>% 
  filter(power >= 0.8) %>% 
  slice_min(n)

# Párový t-test ------------------------------------------------------

# Spočítejte, jako velikost vzorku bychom potřebovali pro dosažení 
# aspoň 80% statistické síly u oboustranného párového t-testu s hladinou 
# alfa = 0,05, ale pro různé velikostí účinku Cohenovo d = 0,1 až 0,8 
# (po krocích o velikosti 0,01). Vhodně výsledky vizualizujte.

pwr_paired_ttest <- pwr_analysis(
  list(d = seq(0.1, 0.8, by = 0.01)),
  wp.t,
  power = 0.8,
  type = "paired")

pwr_paired_ttest

pwr_plot(pwr_paired_ttest,
         x = "d",
         y = "n")

pwr_paired_ttest <- pwr_analysis(
  list(d = seq(0.1, 0.8, by = 0.01),
       power = seq(0.7, 0.9, by = 0.1)),
  wp.t,
  type = "paired")

pwr_paired_ttest

pwr_plot(pwr_paired_ttest,
         x = "d",
         y = "n",
         color = "power")

# Nezávislý t-test ------------------------------------

# Spočítejte, jaké statistické síly pro oboustranný nezávislý t-test s 
# hladinou alfa = 0,05 bychom dosáhli při velikosti vzorku od 5 do 500 osob 
# v každé skupině (předpokládejte početně vyrovnané skupiny) a pro velikosti 
# účinku Cohenovo d = 0,1, 0,2, 0,3 atd. až 0,8. Výsledky vhodně vizualizujte.
pwr_ind_ttest <- pwr_analysis(
  list(n1 = 5:500,
       d = seq(0.1, 0.8, by = .1)),
  wp.t, 
  type = "two.sample"
)


pwr_ind_ttest

pwr_plot(pwr_ind_ttest,
         x = "n",
         y = "power",
         color = "d",
         hline = 0.8)

# Paramtest --------------------------------------------
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
                  alpha = .05) { # arugment pro stanovení hladiny alfa
  
  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 <- t$statistic # uložení testové statistiky t
  p <- t$p.value # uložení p-hodnoty
  
  return(c(t = stat, p = p, sig = (p < alpha))) # 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 <- 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 = 6, # počet jader procesoru k využití
  var.equal = TRUE,
  m1 = 0, m2 = 0.5, sd1 = 0.5, sd2 = 1.5) # ostatní vstupní argumenty

# Zkuste generovat data za předpokladu, že nulová hypotéza platí 
# (obě populace mají shodný průměr rovný např. m1 = m2 = 0). 
# První populce má ale mnohem vyšší směrodatnou odchylku (např. sd1 = 5) 
# než druhá (sd2 = 1). Zkuste, jestli bude očekávaný podíl signifikantních 
# výsledků odpovídat pozorovanému (simulovanému), i když použijeme t-test 
# předpokládající shodu rozptylů (var.equal = TRUE). Protože nulová 
# hypotéza platí, očekávali bychom, že pouze 5% p-hodnot bude menších než 0,05. 
# Zkuste to ověřit pro různý poměr velikostí skupin n1/n2: 10/90, 25/75 a 50/50.
# Pak zkuste použít Welchův t-test, nepředpokládající shodu rozptylů

pwr_var_equal <- grid_search(
  t_fun, # použitá funkce
  params = list(n1 = c(10, 25, 50), n2 = c(90, 75, 50)), # 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 = 6, # počet jader procesoru k využití
  var.equal = TRUE, # shoda rozptylů
  m1 = 0, m2 = 0, sd1 = 5, sd2 = 1)

pwr_var_equal <- results(pwr_var_equal) %>% 
  as_tibble()

pwr_var_equal

pwr_var_equal %>% 
  group_by(n1.test, n2.test) %>% 
  summarise(p_sig = mean(sig))

pwr_var_unequal <- grid_search(
  t_fun, # použitá funkce
  params = list(n1 = c(10, 25, 50), n2 = c(90, 75, 50)), # 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 = 6, # počet jader procesoru k využití
  var.equal = FALSE, # 
  m1 = 0, m2 = 0, sd1 = 5, sd2 = 1) # ostatní vstupní argumenty


pwr_var_unequal <- results(pwr_var_unequal)

pwr_var_unequal %>% 
  group_by(n1.test, n2.test) %>% 
  summarise(p_sig = mean(sig))

# Simulace na základě korelační matice ---------------------------------

lm_fun <- function(simNum, # musí obsahovat argument simNum pro stanovení počtu simulací
                   n, # velikost vzorku
                   cov_mat, # kovariační nebo korelační matice 
                   mean_vector = NULL,# Vektor s průměry
                   formula) { # rovnice lineární regrese 
  
  # Pokud je mean vektor prázdný, předpokládat standardizované proměnné
  # s průměrem nula
  if (is.null(mean_vector)) {
    mu <- rep(0, ncol(cov_mat))
  } else {
    mu <- mean_vector
  }
  
  # Nejdříve si generujeme data z multivariačně normálního rozdělení
  data <- MASS::mvrnorm(n = n, mu = mu, Sigma = cov_mat)
  
  # Fitujeme lineárně regresní model
  fit <- lm(formula, data = as.data.frame(data))
  
  # Extrahujeme koeficienty
  # Funkce k extrakci koeficientů do jednoho atomického vektoru
  get_coefs <- function(fit) {
    out <- fit |> 
      summary() |> 
      coef() |> 
      as.data.frame() |> 
      setNames(c("b", "se", "t", "p"))
    
    row_names <- rownames(out)
    row_names[row_names == "(Intercept)"] <- "int"
    
    out <- out |> 
      dplyr::mutate(predictor = row_names) |> 
      tidyr::pivot_wider(names_from = predictor, values_from = -predictor)
    
    nms <- colnames(out)
    
    out <- as.numeric(out)
    names(out) <- nms
    
    out
  }
  
  # Výstupní hodnoty
  return(get_coefs(fit))
}

# Chceme pomocí lineární regrese zkoumat efekt veku (vekr), 
# školního prospěchu (gpa), monitorování rodiči (monit), 
# důvěry vůči vrstevníkům (duv_v) a optimismu na životní spokojenost (ziv_sp). 
# Předpokládáme, že tyto proměnné mají přibližně multivariačně normální rozdělení 
# a že variance a kovariance mezi tětito proměnnými jsou takovéto:

# Variačně-kovariační matice
cov_mat <- structure(c(0.23, 0.02,-0.07, 0.10, 0.00, 0.81, 
                       0.02, 3.96, 0.56,-0.30, 0.18,-0.01, 
                       -0.07, 0.56, 0.56,-0.10, 0.01,-0.14, 
                       0.10,-0.30,-0.10, 0.32,-0.01, 0.30, 
                       0.00, 0.18, 0.01,-0.01, 0.26, 0.07, 
                       0.81,-0.01,-0.14, 0.30, 0.07, 9.71),
                     dim = c(6L, 6L), 
                     dimnames = list(c("ziv_sp", "vekr", "gpa", "monit", 
                                       "duv_v", "optim"), 
                                     c("ziv_sp", "vekr", "gpa", "monit", 
                                       "duv_v", "optim")))
cov_mat

# Vektor s průměry
mean_vec <- c(ziv_sp = 2.9, vekr = 14.11, gpa = 2.2, monit = 2.8, 
              duv_v = 2.7, optim = 16.32)
mean_vec


# Chceme zjistit statistickou sílu pro testy egresních koeficientů jednotlivých
# prediktorů při velikost vzorku 30, 100, 300, 500 nebo 1000 osob.

pwr_lm <- grid_search(
  lm_fun, 
  params = list(n = c(30, 100, 300, 1000)), # Velikosti vzorku
  n.iter = 2000, # počet simulací
  output = 'data.frame', # Formát výstupu
  cov_mat = cov_mat, # kovariační matice
  mean_vector = mean_vec,
  formula = ziv_sp ~ vekr + gpa + monit + duv_v + optim, # regresní rovnice
  parallel = 'snow', # Paralelní výpočty 
  ncpus = 6 # POčet použitých jader procesoru
) 

pwr_lm <- results(pwr_lm) %>% 
  as_tibble()

pwr_lm

# Sloupce p-hodnot
pwr_lm_pvals <- pwr_lm %>% 
  select(iter, n.test, starts_with("p_"))

# Převod na logner formát
pwr_lm_pvals <- pwr_lm_pvals %>% 
  pivot_longer(-c(iter, n.test), 
               names_to = "predictor", 
               values_to = "p_value") 

# Odhad statistické síly (což je vlastně podíl signifikantních výsledků)
pwr_lm_pvals <- pwr_lm_pvals %>% 
  mutate(sig = p_value < 0.05) %>% # hladina alfa
  group_by(n.test, predictor) %>% 
  summarise(p_sig = mean(sig))

pwr_lm_pvals

# Power curves
pwr_lm_pvals %>% 
  pwr_plot(x = "n.test",
           y = "p_sig",
           color = "predictor")