library(rstatix)
library(effectsize)
library(rcompanion)
library(tidyverse)

# Import dat ------------------------------------------------------------
df <- haven::read_sav("data/eukids.sav")

# Převedení některých proměnných na faktory
df <- df %>% 
  mutate(
    sex = factor(pohlavi,
                 levels = c(1, 2),
                 labels = c("Boys", "Girls")),
    cohort = factor(kohorta,
                    levels = c(6, 10),
                    labels = c("Younger", "Older")),
    .after = id
  )

# Pohlaví - četnosti
df %>% 
  count(sex) %>% 
  mutate(p = n/sum(n))

df %>% 
  count(cohort) %>% 
  mutate(p = n/sum(n))

# Odstranění respondentů s chybějícími hodnotami v sex a cohort
df %>% 
  count(sex) # Tři missing values

df %>% 
  count(cohort) # Jen jedna missing value

df <- df %>% 
  drop_na(sex, cohort)

# Jenom přemístění sloupců někde na začátek
df <- df %>% 
  relocate(monit, warm_m, warm_o,
           .after = cohort)

df
# Srovnání dvou nezávislých výběrů ------------------------------
# Budeme analyzovat rozdíly mezi kohortami v monitorování rodiči

# Chybějících dat moc není
df %>% 
  count(missing = is.na(monit))

# * Deskriptivní statistiky ------------------------------------

# funkce pro výpočet deskriptivních statistik
my_stats <- function(x, level = .95) {
  x <- x[!is.na(x)]
  n <- length(x)
  m <- mean(x)
  s <- sd(x)
  se <- s/sqrt(n)
  p <- level + (1-level)/2
  me <- qt(p, df = n - 1)*se
  
  q <- quantile(x, probs = c(0, .25, .50, .75, 1),
                na.rm = TRUE)

  tibble(m = m, # průměr
         m_lo = m - me, # spodní limit CI
         m_hi = m + me, # horní limit CI
         n = n, # počet validních hodnot
         se = se, # standardní chyba
         skew = psych::skew(x), # šikmost
         kurt = psych::kurtosi(x), # špičatost
         min = q[1], # minimum
         q25 = q[2], # první kvartil
         mdn = q[3], # medián
         q75 = q[4], # třetí kvartil
         max = q[5]) # maximum
}

my_stats(df$monit)

# Deskriptivní statistiky pro monitorování podle kohorty
monit <- df %>% 
  group_by(cohort) %>% 
  summarise(my_stats(monit))

monit

# * Testy ověřujíc některé předpoklady --------------------------

# Levenův test shody rozptylů
df %>% 
  rstatix::levene_test(monit ~ cohort)

# Shapirův-Wilkův test of normality
df %>% 
  group_by(cohort) %>% 
  rstatix::shapiro_test(monit)

# * Grafy ----------------------------------------------

# Jenom změna estetického schématu
theme_set(
  theme_classic(base_size = 15)
)

# Boxplot
df %>% 
  ggplot(aes(cohort, monit)) +
  geom_boxplot()

# Boxplot i s průměry
df %>% 
  ggplot(aes(cohort, monit)) +
  geom_boxplot() +
  stat_summary(fun.data = mean_cl_normal, # funkce pro výpočet CI pro průměr
               geom = "pointrange", # Použitý geom
               color = "red")

# Boxplot + violinplot
df %>% 
  ggplot(aes(cohort, monit)) +
  geom_violin(width = .5) +
  geom_boxplot(width =.25)

# Line charts of means
df %>% 
  ggplot(aes(cohort, monit)) +
  stat_summary(fun.data = mean_cl_normal,
               geom = "pointrange",
               size = 1)

# Kdybychom neměli dopředu vypočtené deskriptivní statistiky
df %>% 
  ggplot(aes(cohort, monit)) +
  stat_summary(fun.data = mean_cl_normal, # vypočíst CI pro průměr
               geom = "errorbar", # Chybové úsečky pro CI průměru
               width = .1) +
  stat_summary(fun = mean, # vypočíst průměr
               geom = "point", # tečkou vyznačit průměr
               size = 5) +
  stat_summary(fun = mean, # vypočíst průměr
               geom = "line", # Spojit body čárou
               mapping = aes(group = 1)) # Jen jednou

# Stejný graf, ale využívající dopředu vypočtených deskriptivních statistik
monit %>% 
  ggplot(aes(x = cohort,
             y = m)) +
  geom_point(size = 5) +
  geom_line(mapping = aes(group = 1)) +
  geom_errorbar(mapping = aes(ymin = m_lo,
                              ymax = m_hi),
                width = .1)

# Jittered points 
df %>% 
  ggplot(aes(cohort, monit)) +
  geom_jitter(width = .1,
              height = .05,
              size = 5,
              alpha = .2)

# Ověření normality

# Převod na z-skóry
scaled <- df %>% 
  select(cohort, monit) %>% 
  group_by(cohort) %>% 
  mutate(
    zscores = as.double(scale(monit))
  )


# Histogramy a grafy hustoty pravděpodobnosti
scaled %>% 
  ggplot(aes(zscores)) +
  geom_histogram(
    mapping = aes(y = after_stat(density)), # Na ose Y hustota pravděpodobnosti místo absolutních četností
    color = "black", fill = "grey"
  ) +
  geom_density(linewidth = 1, 
               color = "blue") +
  stat_function(fun = dnorm, # Zakrelit i křivku normálního rozdělení
                color = "red",
                linewidth = 1,
                linetype = "11") +
  facet_wrap(~cohort, nrow = 2)

# Q-Q grafy
scaled %>% 
  ggplot(aes(sample = zscores)) +
  geom_qq_line() +
  geom_qq() +
  facet_wrap(~cohort)

# * Statistické testy pro srovnání dvou nezávislých výběrů -----------------

# Nezávislý t-test 
t.test(monit ~ cohort, 
       data = df)

t.test(monit ~ cohort, 
       data = df,
       var.equal = TRUE) # Předpokládat shodný rozptyl

t.test(monit ~ cohort, 
       data = df,
       conf.level = .99, # Nastavení úrovně spolehlivosti pro CI
       var.equal = TRUE) # Předpokládat shodný rozptyl

# Wilcox-Mann-Whitney test (neparametrická alternativa nezávislého t-testu)
wilcox.test(monit ~ cohort, 
            data = df,
            conf.int = TRUE)

# * Velikosti účinku -------------------------------------------

# Cohenovo d
effectsize::cohens_d(monit ~ cohort, 
                     data = df)
effectsize::cohens_d(monit ~ cohort, 
                     pooled_sd = FALSE, # Předpoládat shodu rozptylů?
                     ci = .99, # Nastavení úrovně spolehlivosti
                     data = df)


# Velikost účinku: Pravděpodobnost superiority
rcompanion::wilcoxonPS(monit ~ cohort,
                       ci = TRUE,
                       data = drop_na(df, monit))

# Srovnání dvou závislých výběrů --------------------------------
# Budeme testovat rozdíl mezi vnímanou vřelostí matky a otce

df %>% 
  count(warm_o_missing = is.na(warm_o)) # Více chybějích hodnot než u warm_m

df %>% 
  count(warm_m_missing = is.na(warm_m))

# Příprava dat
df_clean <- df %>% 
  drop_na(warm_m, warm_o) %>% 
  select(id, warm_m, warm_o)
df_clean
# Převod na delší formát
df_long <- df_clean %>% 
  pivot_longer(warm_m:warm_o,
               names_to = "variable")
df_long

# * Deskriptivní statistiky a test normality rozdílových skórů -----------
warmth <- df_long %>% 
  group_by(variable) %>% 
  summarise(my_stats(value))
warmth


df_clean %>% 
  mutate(diff = warm_m - warm_o) %>% 
  rstatix::shapiro_test(diff)

# Grafy ---------------------------------------------------------------

# Boxplot 
df_long %>% 
  ggplot(aes(variable, value)) +
  geom_boxplot()

# Boxplot + violin plot
df_long %>% 
  ggplot(aes(variable, value)) +
  geom_violin(width = .50) +
  geom_boxplot(width = .25)

# Line graphs of means
warmth %>% 
  ggplot(aes(x = variable,
             y = m,
             ymin = m_lo,
             ymax = m_hi,
             group = 1)) +
  geom_point(size = 5) +
  geom_line() +
  geom_errorbar(width = .1) +
  scale_y_continuous(breaks = seq(3, 4, by = .05))

# Bodový graf 1:1
df_clean %>% 
  ggplot(aes(warm_m, warm_o)) +
  geom_jitter(height = .03,
              width = .03,
              alpha = .25,
              size = 3) +
  coord_equal(xlim = c(1, 4),
              ylim = c(1, 4)) +
  geom_abline()

# Histogram a graf hustoty pravděpodobnosti pro rozdílové skóry   
df_clean <- df_clean %>% 
  mutate(diff = warm_m - warm_o)

df_clean %>% 
  ggplot(aes(x = diff)) +
  geom_histogram(mapping = aes(y = after_stat(density)),
                 fill = "grey",
                 color = "black") +
  geom_density(linewidth = 1, 
               color = "blue") +
  stat_function(fun = dnorm, # Zakreslit i křivku normálního rozdělení
                color = "red",
                linewidth = 1,
                linetype = "11",
                args = list(mean = mean(df_clean$diff), # nastavit stejný průměr a směrodatnou odchylku jako data
                            sd = sd(df_clean$diff)))

# Q-Q graf
df_clean %>% 
  ggplot(aes(sample = diff)) +
  geom_qq_line() +
  geom_qq()

# Těžké chvosty jsou jasně vidět z obou předchozích grafů

# * Statistické testy pro srovnání dvou závislých výběrů -----------------

# Párový t-test 
t.test(df_clean$warm_m, df_clean$warm_o,
       paired = TRUE)

t.test(df_clean$warm_m, df_clean$warm_o,
       paired = TRUE,
       data = df_clean,
       conf.level = .99)

# Alternativní zápis
t.test(Pair(warm_m, warm_o) ~ 1,
       data = df_clean)


# Wilcoxonův test (neparametrická alternativa závislého t-testu)
# Pracuje s pořadím (ranks)
wilcox.test(df_clean$warm_m, df_clean$warm_o,
            paired = TRUE)

wilcox.test(Pair(warm_m, warm_o) ~ 1,
            data = df_clean)


# Znaménkový test: sleduje pouze to, kolikrát  warm_m > warm_o nebo warm_m < warm_o
df_clean %>% 
  count(warm_m_higher = warm_m > warm_o)

df_clean %>% 
  count(nontied = warm_m != warm_o)

binom.test(x = 373, # kolikrát warm_m > warm_o
           n = 576) # celkový počet non-ties


# Velikost účinku --------------------------------
# Cohenovo d
effectsize::rm_d(df_clean$warm_m, df_clean$warm_o,
                 paired = TRUE,
                 ci = .95)

# Pravděpodobnost superiority
# Spočítal už vlastně binomický test (jako probability of success)
binom.test(x = 373, # kolikrát warm_m > warm_o
           n = 576) # celkový počet non-ties