##%######################################################%##
#                                                          #
####                 Statistical Power                  ####
#                                                          #
##%######################################################%##


# Packages ----------------------------------------------------------------
library(WebPower)
library(lme4)
library(paramtest)
library(simstandard)
library(lavaan)
library(tidyverse)
library(mixedpower)
library(simr)

# WebPower -------------------------------------------------------------
# balíček WebPower poskytuje je takovou eRkovou obdobou GPower
# v tom smyslu, že umožňuje základní odhad statistické síly pro 
# základní statistické testy

# * Korelace ------------------------------------------------------------
# Nejdříve si vyzkoušíme odhad statistické síly pro korelace
# Dejme tomu, že zkoumáme vztah mezi stresem a well-beingem
# Na základě předchozích výzkumů očekáváme korelaci r = 0.3.
# Můžeme si položit otázku: kdybychom získali měřili stres a well-being
# u 50 osob, jakou statistickou sílu bychom měli při alfa = 0.05

# můžeme 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)

# defaultně se počítá s oboustranným testem
# ale můžeme použít argument alternative pro stanovení "směru" alternativní
# hypotézy
# Také můžeme změnit nulovou hypotézu pomocí argumentu rho0 (defalt je rho0 = 0)

# Kdybychom chtěli vypočíst statistckou sílu pro různé velikosti vzorku
wp.correlation(n = c(10, 20, 30), # Jeden z argumentů může zahrnovat více hodnot
               r = 0.3, alpha = .05)

# Když uvedeme power, ale nikoli n, vypočte se potřebná velikost vzorku
wp.correlation(power = .8, r = 0.3, alpha = .05)

# Partial correlation
wp.correlation(power = .8, r = 0.3, alpha = .05, 
               p = 1) # počet "kontrolovaných" proměnných

# Co když chceme zjistit power pro různé velikosti vzorku, různé r a různou alfa
# Vytvoříme si list() relevantních hodnot
df_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
)

# Funkcí cross_df() získáme všechny potřebné kombinace jako řádky dataframu
df_crossed <- cross_df(df_input)
df_crossed

# Pomocí funkce pwr_df pak dostaneme výsledky ve formě dataframu
# argumenty jsou data a funkce z balíčku WebPower
df_final <- pwr_df(df_crossed, fn = wp.correlation)
df_final

# Výsledný dataframe pak můžeme použít pro tvorbu grafu
df_final %>% 
  ggplot(aes(n, power, color = as.factor(r))) +
  geom_line(size = 1) +
  facet_wrap(~alpha, labeller = "label_both") +
  geom_hline(yintercept = .8, linetype = 2) +
  theme(
    legend.position = "top"
  )

# Taky samozřejmě můžeme chtít zjistit potřebnou velikost vzorku pro
# stanovenou velikost účinku, sílu a alfa
df_input <- list(
  r = seq(0.2, 0.6, by = 0.01),
  alpha = c(0.05, 0.01),
  power = c(.8, .9, .95)
)

df_final <- df_input %>% 
  cross_df() %>% 
  pwr_df(fn = wp.correlation)
df_final

df_final %>% 
  ggplot(aes(r, n, color = as.factor(power))) +
  geom_line(size = 1) +
  facet_wrap(~alpha) +
  labs(x = "Pearson r", y = "Sample size", color = "Power") +
  theme(legend.position = "top") +
  geom_vline(xintercept = .3, linetype = 2)
  
# * t-tests -----------------------------------------------------------------

### Párový t-test
# Dejme tomu, že chceme testovat efekt nějaké intervence; máme skupinu
# studentů a měříme je před a po intervenci (čili párová data)
# Očekáváme velikost účinku d = 0.3
# Kolik účastníků potřebujeme k dosažení statistické síly 0.8
wp.t(d = .3, power = 0.8, alpha = 0.05, type = 'paired')

### Nezávislý t-test
# Dejme tomu, že místo vnitrosubjektového použijeme mezisubjektový desing,
# takže máme dvě nezávislé skupiny, jedna bude vystavena placebo-intervenci
# a druhá skutečné intervenci
# Očekáváme stejnou velikost účinku a chceme mít vyvážený počet osob
# v obou skupinách
wp.t(d = .3, power = 0.8, alpha = 0.05, type = 'two.sample')


# Co kdyby ale poměr skupin nebyl vyvážený
# Např. kdybychom chtěli poskytnout skutečnou intervenci co nejvíce lidem,
# ale věděli bychom, že celková velikost vzorku, kterou můžeme získat, 
# je 500 osob. Stále bychom ale chtěli mít sílu testu aspoň 0.8
# Jak moc mohou být skupiny nevyvážené?

df_input <- tibble(
  n1 = 2:498,
  n2 = 500 - n1,
  d = 0.3,
  alpha = 0.05,
  type = 'two.sample.2n'
)

df_output <- pwr_df(df_input, fn = wp.t)
df_output

df_output %>% 
  ggplot(aes(n1, power)) +
  geom_line() +
  geom_hline(aes(yintercept = .8)) +
  labs(x = "n1 (n2 = 500 - n1)")

df_output %>% 
  filter(power >= .8) %>% 
  slice_min(n1)


# * ANOVA -------------------------------------------------------------------

### One-Way ANOVA
# Dejme tomu, že chceme srovnat čtyři skupiny a očekáváme velikost účinku
# R2 = 0.2
# Kdybychom mohli získat pouze pro každou skupinu pouze 10 účastníků
# (tj. dohromady 40 osob), jakou statistickou sílu by měl F-test ověřující
# shodu průměrů skupin.

# Nejdříve vypočteme Cohenovo f: 
# odmocninu z (přírůstku) vysvětleného rozptylu dělená podílem
# nevysvětleného rozptylu po zahrnutí daného efektu
#  sqrt((r2_full - r2_null) / (1 - r2_full))

cohen_f <- r2_to_f(.2)

wp.anova(k = 4, n = 40, alpha = .05, f = cohen_f)

# Také můžeme vypočíst minimální velikost účinku pro požadovanou 
# statistickou sílu, velikost vzorku a hladinu alfa

# One can also calculate the minimum detectable effect to achieve certain power 
wp.anova(f = NULL, k = 4, n = 100, power = 0.8, alpha = 0.05)$f %>% 
  f_to_r2()

# Anebo požadovanou velikost vzorku pro danou velikost účinku, alfu a počet 
# skupin
wp.anova(f = cohen_f, k = 4, n = NULL, power = 0.8, alpha = 0.05)

# Funkce wp.kanova je zobecněním pro více nezávislých proměnných (k-way ANOVA)
wp.kanova(n = NULL,  
          power = .8, # požadovaná statistická síla
          ndf = 3, # stupně volnosti pro daný efekt, pro hlavní efekty: počet skupina - 1
          ng = 4, # celkový počet skupin (násek počtu úrovní jednotlivých faktorů)
          f = cohen_f, # velikost účinku
          alpha = .05)

# Převod z parciální eta2 či omega2 na Cohenovo f:
# sqrt(peta2 / (1 - peta2))

peta2_to_f(0.3)

# * Linear regression -------------------------------------------------------

# Výzkumník se domnívá, že průměr známek ze střední školy a výsledek TSP
# vysvětluje 25 % rozptylu v průměru známek na vysoké škole.
# Pokud má vzorek 50 osob, jaká je velikost účinku pro tento efekt
# (alha = 0.01)
f2 <- r2_to_f(0.25, f_squared = TRUE)
f2

wp.regression(n = 50, 
              p1 = 2, # počet preditkorů
              f2 = f2,
              alpha = 0.01)

# Jiný výzkumník se domnívá, že hodnocení z přijímacího pohovoru
# vysvětluje dalších 5 % rozptylu.
# Jakou velikost vzorku potřebujeme získat pro dosažení 80% statistické 
# síly při alpha = 0.01
f2 <- r2_to_f(0.3, r2_null = 0.25, f_squared = TRUE)
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
              alpha = 0.01) # Hladina alfa


# simr and mixedpower packages -------------------------------------------------

# * Scénář 1: Použití stávajících dat -------------------------------------
# Yan et all. (2014) tested 48 subjects, 
# each of whom read 120 sentences investigating the effect of different factors 
# on the first landing position (FLP, i.e. the position in a sentence 
# your eyes first land on).
load("data/Yan_et_al.RData") 
glimpse(YanData)


# Autoři hypotetizovali efekt délky slov a složitosti věty
# Odhadujeme také náhodný průsečík napríč respondenty a větami
# Protože lze očekávat, že se průměrná FLP se liší člověk od člověka
# a také v závislosti na konkrétní větě i po kontrole délky a složitosti
flp_model <- lmer(flp ~ word_length * complexity + (1|subject) + (1|sentence),
                  data = YanData)

summary(flp_model, corr = F) # let's have a look!

# ** Různé velikosti vzorku -----------------------------------------------
# Protože chceme zjistit statistickou sílu pro různé velikosti vzorku,
# což odpovídá náhodné proměnné "subject",
# můžeme použít funkci mixedpower() s následujícími argumenty

pwr_flp <- mixedpower(
  model = flp_model, # použitý mixed model
  data = YanData, # Dataset pro fitování modelu
  fixed_effects = c("word_length", "complexity"), # Všechny fixní efekty
  simvar = "subject", # Náhodný efekt, který chceme variovat
  steps = c(30, 60), # Ú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
  databased = TRUE) # Odhad statistické síly na základě efektů vypočtených z dat?

multiplotPower(pwr_flp)

# ** Různé velikosti vzorku a účinku ------------------------------------------
# Pomocí arugmentu sesoi můžeme specifikovat vlastní velikosti účinku
# (jiné, než byly odhadnuty na základě dat)
# protože můžeme chtít odhad statistické síly i pro menší, ale stále 
# prakticky zajímavé efekty, a protože např. kvůli publikačnímu zkreslení
# bývají publikované efekty nadhodnocené

summary(flp_model)
# Co kdybychom chtěli zjistit statistickou sílu i pro menší efekty?
flp_model
# Specifikace smallest effect sizes of interest
sesoi <- c(3.66, 0.75, -0.065, 0.09) 
# v pořadí dle modelu
flp_model
# První je intercept, pak word_length, comlexity a nakonec jejich interakce
pwr_subjects <- mixedpower(
  model = flp_model, # použitý mixed model
  data = YanData, # Dataset pro fitování modelu
  fixed_effects = c("word_length", "complexity"), # Všechny fixní efekty
  simvar = "subject", # 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
  databased = TRUE, # Odhad statistické síly (i) na základě efektů z dat?
  SESOI = sesoi # Odhad statistické síly pro stanovené velikosti účinku
  )
dev.off()
multiplotPower(pwr_subjects, filename = "pwr.png")

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

# ** Různý počet stimulů --------------------------------------------------
# Můžeme samozřejmě měnit nejen velikost vzorku, ale také počet stimulů
# (v tomto případě počet vět, které každý účasnítk četl)
pwr_sentences <- mixedpower(
  model = flp_model,
  data = YanData, 
  fixed_effects = c("word_length", "complexity"),
  simvar = "sentence",
  steps = c(100, 150, 200),
  critical_value = 2,
  n_sim = 100,
  databased = TRUE
)

multiplotPower(pwr_sentences, filename = "plots/pwr_sentences.png")

# ** Různy počet participantů a stimulů -------------------------------------
# Pomocí funkce R2power můžeme variovat dva náhodné efekty
# ale v současnosti funkce vyžaduje, aby argument R2level měl pouze jednu 
# hodnotu, např. jednu konkrétní velikost vzorku

pwr_R2 <- R2power(model = flp_model, 
                  data = YanData,
                  fixed_effects = c("word_length", "complexity"),
                  simvar = "sentence", 
                  steps = c(100, 150, 200),
                  R2var = "subject", 
                  R2level = 30, # v současné době podpuruje funkce pouze jednu hodnotu druhého náhodného efektu
                  critical_value = 2,
                  n_sim = 50,
                  SESOI = sesoi,
                  databased = T)


# * Scénář 2: Simulace dat ---------------------------------------
# Musíme si tedy nejdříve vytvořit dummy data
# Dejme tomu, že chceme testovat žáky ve třech časových bodech.

# Máme přístup k 5 školním třídám a plánujeme z každé třídy vybrat 10 žáků
variables <- list(
  pupil = 1:10,
  class = 1:5,
  time = 0:2
)
df <- cross_df(variables) # vytvoření všech kombinací
df %>% View()

# Polovinu žáků z každé třídy chceme vystavit experimentální 
# a polovinu kontrolní podmínce
df <- df %>% 
  mutate(
    group = c("exp", "ctrl") %>% 
      rep(length.out = nrow(df)) %>% 
      factor()
  )
df %>% View()

# V dalším kroku musíme specifikovat parametry modelu
# Použijeme tento model
# y ∼ group + time + group*time + (1 | class / id) + residual
# který zahrnuje efekt skupiny, času a jejich interakce
# plus náhodný průsečík pro žáky "zanořené" do tříd
# Náhodné průsečík pro třídy modeluje to, že se jednotlivé třídy mohou
# lišit v průměrné úrovni závislé proměnné
# Náhodný průsečík pro žáky modeluje to, že se liší i jednotliví žáci 
# v rámci jedné třídy v průměrné úrovni závislé proměnné 

# Pro demonstrační účely použijeme arbitrární hodnoty,
# Ale kdybychom plánovali skutečný výzkum, museli bychom je pečlivě zvolit

# Nejdříve stanovíme fixní efekty
fixed_ef <- c(
  5, # intercept 
  0, # efekt intervence na začátků výzkumu (kdy vlastně ještě žádná intervence neproběhla)
  0.2, # efekt času v kontrolní skupině
  0.7 # efekt interakce, o kolik více porostou skóry ve skupině s intervencí
)

# Pak náhodné efekty
rand_ef <- list(
  0.5, # rozptyl vysvětlený rozdíly mezi žáky v rámci třídy
  0.1) # rozptyl vysvětlený rozdíly mezi třídami

# Residual variance
res <- 2 # reziduální rozptyl

# Pomocí funkce makeLmer() z balíčku simr pak zkombinujeme všechny parametry
# modelu dohromady
simulated_model <- makeLmer(y ~ group*time + (1 | class / pupil), 
                            fixef = fixed_ef, VarCorr = rand_ef, 
                            sigma = res, data = df)
summary(simulated_model) # Vidíme, že všechny efekty jsou takové,
                         # jako jsme si stanovili 

# Následně můžeme použít funkci mixedpower() nebo R2power()
# pro odhad statistické síly při různém počtu žáků nebo tříd,
# případně můžeme pomocí sesoi stanivit jiné velikosti účinku
simulated_pwr <- mixedpower(model = simulated_model, 
                            data = df,
                            fixed_effects = c("group", "time"),
                            simvar = "class", 
                            steps = c(5, 10, 20),
                            critical_value = 2, 
                            n_sim = 100)


# Balíček paramtest -------------------------------------------------------
# Nejflexibilnější přístup pro power analýzu pomocí simulace,
# Ale zároveň nejsložitější pro naprogramování

# Začněme jednoduše výpočtem power pro nezávislý t-test
# Samozřejmě bychom mohli použít přesné, analytické řešení bez simulace
# např. pomocí wp.t() funkce z WebPower
wp.t(n1 = 20, n2 = 40, d = 0.5, alpha = .05, type = "two.sample.2n")

# Ale tento přístup má omezení, protože např. předpokládá shodu rozptylů skupin

# * Nezávislý t-test a funkce run_test() ------------------------------------

# paramtest balíček nejprve vyžaduje napsat si vlastní funkci pro generování
# dat, provedení daného statistického testu a uložení výsledků

# První argument by měl být vžy simNum, ten slouží ke stanovení počtu simulací
# další argumenty pak upřesňují, jaká data budeme generovat
# a závisejí na tom, které funkce k tomu použijeme a jaký test chceme použít

# Obecně musí funkce mít následující strukturu
  # 1. Kód sloužící ke generování dat
  # 2. Kód sloužící k odhadu (fitu) modelu/ uložení objektu s fitem
  # 3. Kók sloužící k extrakci potřebných informací z fit-objektu
    # obvykle se jedná o klíčové koeficienty společně s p-hodnotami apod.
  # 4. Stanovení výstupních hodnot v podobě atomického numerického vektoru
    # jehož prvky mají jména

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
  p <- t$p.value # uložení p-hodnoty
  
  return(c(t = stat, p = p, sig = (p < alpha))) # 
  # 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
}

# Pro provedení simulace samotné pak použijeme funkci run_test()
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
  n1 = 20, n2 = 40, m1 = 0, m2 = 0.5, sd1 = 0.5, sd2 = 1.5) 
  # nakonec argumenty doplněné do funkce t_func
  # chceme simulovat data pro dvě skupiny o počtu 20 a 40 osob
  # předpokládáme populační průměry 0 a 0.5
  # a populační sd 0.5 a 1.5
  # arugmenty var.equal a alpha mají defaultní hodnoty,
  # můžeme je ale klidně změnit, kdybychom chtěli použít klasický
  # t-test s rovností rozptylů nebo jinou hladinu alfa

# Výsledky extrahujeme pomocí funkce results()
# každý řádek je jedna simulace
results(pwr_ttest) %>% head()

# Odhad power pak vypočteme jednoduše jako podíl signifikantních testů
results(pwr_ttest) %>%
  summarise(power = mean(sig),
            sum_sig = sum(sig))

# Kdybychom chtěli vypočíst interval spolehlivosti,
# můžeme použít funkci binom.test()
binom.test(x = 944, # počet "úspěchů" (signifikantních výsledků) 
           n = 2000) # počet pokusů (celkový počet simulací)


# * Různé vstupní hodnoty a funkce grid_search() --------------------------

# V předchozím případě jsem použili funkci run_test(), protože jsme
# jsme použili jenom jedny vstupní hodnoty.
# Pomocí funkce grid_search jich ale můžeme použít několik 

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í
  m1 = 0, m2 = 0.5, sd1 = 0.5, sd2 = 1.5) # ostatní vstupní argumenty

res <- results(pwr_ttest) %>%
  group_by(n1.test, n2.test) %>%
  summarise(power = mean(sig))
res

# Power plot pro různé velikosti skupin
res %>%
  ggplot(aes(n1.test, power, color = factor(n2.test) )) +
  geom_point(size = 2) +
  geom_line() +
  labs(x = "n1", color = "n2") +
  coord_cartesian(ylim = c(0,1))



# * Generování dat na základě korelační matice ----------------------------
# Data můžeme simulovat i na základě očekávané populační matice korelací
# Samozřejmě obvykle k jejímu odhadu použijeme vzorek z předchozí studie
# (nejlépe dotatečně velký a reprezentativní)

# Předpokládaná korelační matice
cormat <- matrix(
  c(1.0, 0.3, 0.6,
    0.3, 1.0, 0.4,
    0.6, 0.4, 1.0),
  byrow = TRUE, nrow = 3,
  dimnames = list(
    c("x", "y", "z"),
    c("x", "y", "z"))
)

# Generování multivariačně normálních dat na základě této matice
# Pro zjednodušení pracujeme se standardizovanými proměnnými
# se sigma = 1 a mí = 0.
# Funkce mvrnorm generuje předpokládá multivariačně normální rozdělení
# (z toho generuje data), které je definováno variačně-kovariační maticí

# můžeme si vygenerovat např. 100 hodnot
df <- MASS::mvrnorm(100, 
                    mu = rep(0, 3), # Průměry proměnných
                    Sigma = cormat) # Variačně-kovariační matice


fit <- lm(y ~ x + z, data = as.data.frame(df))
coef(summary(fit))

# Kdybychom chtěli simulovat power pro prediktory v lineárně regresním modelu,
# museli bychom si opět vytvořit vlastní funkci

lm_test <- function(simNum, # musí obsahovat argument simNum pro stanovení počtu simulací
                    n, # velikost vzorku
                    rmat, # korelační matice předpokládaných vztahů
                    formula, # rovnice lineární regrese
                    alpha = 0.05) { # hladina alfa
  
  # Nejdříve si generujeme data z multivariačně normálního rozdělení
  data <- MASS::mvrnorm(n = n, mu = rep(0, nrow(rmat)), Sigma = rmat)
  
  # Fitujeme lineárně regresní model
  fit <- lm(formula, data = as.data.frame(data))
  
  # Extrahujeme koeficienty
  x <- coef(summary(fit))
  x <- x[,-2]
  
  # Převedeme je na atomický vektor se jmény prvků
  y <- as.data.frame(x)
  colnames(y) <- c("b", "se", "p")
  y$sig <- y$p < alpha
  y$pred <- rownames(y)
  y <- tidyr::pivot_longer(y, cols = c(b, se, p, sig), names_to = "what")
  y <- tidyr::unite(y, col = "est", c(pred, what))
  
  out <- y$value
  names(out) <- y$est
  
  # Výstupní hodnoty
  return(out)
}


# A odhadneme power pro jednotlivé prediktory
power_lm <- grid_search(lm_test, 
                        params = list(n = c(50, 100, 200)), 
                        n.iter = 1000, 
                        output = 'data.frame', 
                        rmat = cormat, 
                        formula = y ~ 0 + x + z,
                        parallel = 'snow', 
                        ncpus = 6)

# Shrnutí výsledků
results(power_lm) %>% head()


results(power_lm) %>% 
  ggplot(aes(x_b)) +
  geom_histogram()

results(power_lm) %>%
  group_by(n.test) %>%
  summarise(power_x = mean(x_sig),
            power_z = mean(z_sig))


# * Power pro strukturní modely ---------------------------------------------

# Power analýzu pro celkový test fitu modelu poskytuje balíček WebPower
# Staví na fit-indexu RMSEA
# RMSEA = sqrt(chi2 - df) / sqrt(df*(N-1))
# MacCallum, Browne a Sugawara (1996) původně doporučili cut-offs
# 0.01, 0.05 a 0.08 pro skvělý, dobrý a přijatelný fit

wp.sem.rmsea(n = NULL,
             power = .80, # Požadovaná síla
             df = 12, # Stupně volnosti modelu
             rmsea0 = 0.02, # Nulová hypotéza, nejčastěji jaký misfit na úrovni populace považujeme za zanedbatelný
             rmsea1 = 0.06) # Alternativní hypotéza, nejčastěji nejhorší misfit modelu, který jsme ještě ochotno akceptovat

# Často ale chceme určit power pro klíčové koeficienty modelu
# Za tímto účelem nejdříve musíme specifikovat populační model

pop_m <- "
# Tři faktory (f1 až f3) měřené třemi položkami (celkem devíti, x1 až x9)
f1 =~ 0.8*x1 + 0.9*x2 + 0.7*x3
f2 =~ 0.9*x4 + 0.8*x5 + 0.8*x6
f3 =~ 0.9*x7 + 0.8*x8 + 0.7*x9

# Vztahy mezi faktory
f3 ~ 0.5*f1 + 0.1*f2
f1 ~~ 0.2*f2

# Škála/rozptyl faktorů, nastavena na 1, abychom mohli považovat
# všechny koeficienty za standardizované
f1 ~~ 1*f1
f2 ~~ 1*f2
f3 ~~ 1*f3
"

# Můžeme si vyzkoušet, že funkce generuje data pozorovaných proměnných
my_data <- simulateData(pop_m, sample.nobs = 1000)
head(my_data)

# Pro urychlení výpočtu/rychlejší konvergenci si stanovíme i startovací hodnoty
# Které budou tytéž jako populační hodnoty
start_model <- "
f1 =~ NA*x1 + start(0.8)*x1 + start(0.9)*x2 + start(0.7)*x3
f2 =~ NA*x1 + start(0.9)*x4 + start(0.8)*x5 + start(0.8)*x6
f3 =~ NA*x1 + start(0.9)*x7 + start(0.8)*x8 + start(0.7)*x9

# Parametry, o které máme, označíme libovolnými labely, např. a, b
f3 ~ a*start(0.5)*f1 + b*start(0.1)*f2

# Korelace mezi f1 a f2
f1 ~~ start(0.2)*f2

# Rozptyl faktorů           
f1 ~~ 1*f1
f2 ~~ 1*f2
f3 ~~ 1*f3
"

# Takhle bychom fitovali model na jednom simulovaném datasetu
my_fit <- sem(star_vals, data = my_data)
summary(my_fit, standardized = TRUE)

sem_test <- function(simNum, # počet simulací
                     n, # velikost vzorku
                     pop_model, # populační model (musí obsahovat populační parametry)
                     model, # model se stejnou strukturou, ale bez parametrů (ty použít pouze jako startovací hodnoty)
                     labels, # labely koeficientů, o které se zajímáme
                     alpha = 0.05) { # hladina alfa
  
  # Generovování dat
  data <- lavaan::simulateData(pop_model, sample.nobs = n)
  # Fitování modelu
  fit <- lavaan::cfa(model, data = data)
  # Uložení odhadů parametrů
  ests <- lavaan::parameterEstimates(fit)
  
  # when using parallel processing, we must refer to functions directly
  # Transformace est do podoby atomického vektoru se jmény prvků
  x <- subset(ests, label %in% c("a", "b"))[, c("label", "est", "pvalue")]
  x$sig <- x$pvalue < 0.05
  x <- tidyr::pivot_longer(x, cols = c(est, pvalue, sig))
  x <- tidyr::unite(x, col = "what", name, label)
  
  # Výstupní hodnoty
  out <- x$value
  # Jejich názvy
  names(out) <- x$what
  return(out)
    
}

power_sem <- run_test(
  sem_test,
  params = list(n = c(100, 200, 300)),
  n.iter = 100, 
  output = 'data.frame',
  pop_model = pop_m, 
  model = start_model, 
  labels = c("a", "b"), 
  parallel = 'snow', ncpus = 6)

results(power_sem) %>% 
  group_by(n.test) %>% 
  summarise(
    power_a = mean(sig_a),
    power_b = mean(sig_b)
  )