##%######################################################%##
#                                                          #
####             Základní statistické testy             ####
#                                                          #
##%######################################################%##



# Balíčky -----------------------------------------------------------------
library(tidyverse)
library(ggmosaic)
source("r101_funs.R")

# Příprava dat --------------------------------------------------------

# Import dat
df <- readRDS("data/eukids.rds")

# Přehled dat
glimpse(df)

tibble(
  var_name = names(attributes(df)$variable.labels),
  var_label = attributes(df)$variable.labels
) %>% 
  View()

# Korelace ----------------------------------------------------------------
# Selekce proměnných
variables <- c("optim", "ziv_sp", "selfe", "effi", "zdravi", "deprese")
df_cor <- df %>% select(all_of(variables))
df_cor

# Vytvoření všech možných párů proměnných
combs <- combn(variables, 2) %>% 
  t() %>% 
  as_tibble() %>% 
  setNames(c("x", "y"))

# Histogramy každé proměnné
plot_histogram_mult(df_cor, var = variables)

# Scatterplots
plot_scatter_mult(df_cor,
                  xvar = combs$x, 
                  yvar = combs$y, 
                  geom = "jitter", alpha = 0.25)

# Použití funkce cor z base R
cor(df_cor) # defaultně vyřadí proměnné s chybějícími hodnotami
cor(df_cor, use = "pairwise")

# If use is "everything", NAs will propagate conceptually, 
# i.e., a resulting value will be NA whenever one of 
# its contributing observations is NA.

# If use is "all.obs", then the presence of missing observations will 
# produce an error. 

# If use is "complete.obs" then missing values 
# are handled by casewise deletion (and. 
cor(df_cor, use = "pairwise") %>% 
  round(digits = 2) 
cor(df_cor, use = "pairwise", method = "kendall") %>% 
  round(digits = 2)
cor(df_cor, use = "pairwise", method = "spearman") %>% 
  round(digits = 2)

# funkce corr.test() z balíčku psych
psych::corr.test(df_cor)
psych::corr.test(df_cor, 
                 adjust = "bonferroni", # jakou korekci p-hodnot použít
                 method = "spearman", # typ korelace
                 alpha = 0.005) %>% # hladina spolehlivosti bude 1 - alpha
  print(short = FALSE)

# Custom funkce cor_table()
cor_table(df_cor, 
          adjust = "none", # jakou korekci p-hodnot použít
          conf = .995, # hladina spolehlivosti
          sig_levels = c(0.005, 0.001, 0.0001)) # hvězdičkami značené p-hodnoty

# Korelogramy pomocí balíčku corrplot a stejnojmenné funkce corplot()
rmat <- cor(df_cor, use = "pairwise")
rmat # korelační matice

corrplot::corrplot(rmat)
corrplot::corrplot(rmat, method = 'number', order = "hclust")
corrplot::corrplot.mixed(rmat, 
                         lower = "number",
                         upper = "ellipse",
                         lower.col = "black",
                         order = "hclust")

# Lineární regrese --------------------------------------------------------
# Budeme predikovat depresivitu na základě pohlaví, věku, průměru známek,
# důvěry ve vrstevníky/rodiče a optimismu

# Výběr proměnných a vyčištění dat
df_clean <- df %>% 
  select(id, skola, trida, deprese, pohlavi, vekr, 
         gpa, duv_v, duv_r, selfe, optim) %>% 
  na.omit()

# Průzkum korelací mezi spojitými proměnnými
df_clean %>% 
  select(where(is.numeric), -c(id:trida)) %>% 
  cor_table()

# Vztah kategorického prediktoru a závislé proměnné
df_clean %>% 
  ggplot(aes(pohlavi, deprese)) +
  geom_boxplot()

# Odhad prvního modelu pomocí funkce lm() jako "linear model"
lmfit1 <- lm(deprese ~ pohlavi + vekr, data = df_clean) %>% 
  lm.beta::lm.beta() # doplnění standardizovaných koeficientů pomocí funkce lm.beta() 
            # ze stejnojmenného balíčku

lmfit1 # pouze rovnice a koeficienty
summary(lmfit1) # Celkový přehled modelu

# Odhad druhého modelu
lmfit2 <- lm(deprese ~ pohlavi + vekr + gpa + duv_v + duv_r + selfe + optim, 
             data = df_clean) %>% 
  lm.beta::lm.beta()

# Přírůstek R2
summary(lmfit1)$r.squared
summary(lmfit2)$r.squared
summary(lmfit2)$r.squared - summary(lmfit1)$r.squared

# Srovnání obou modelů pomocí funkce anova()
# nulovou hypotézou je nulový přírůstek vysvětleného rozptylu
anova(lmfit1, lmfit2)

# jaký je rozdíl v R2 mezi modely?
summary(lmfit1) %>% str()

str(lmfit1)
# balíček broom nabízí funkce tidy() a glance() pro zobrazení modelu
# v jedné tibble
  # tidy() extrahuje informace o koeficientech modelu
  # glace() o celkovém fitu modelu

broom::tidy(lmfit2)
broom::glance(lmfit2)

# Pro intervaly spolehlivosti, zde 99.5%
confint(lmfit2, level = 0.995)

# Pro základní grafickou diagnostiku modelu stačí na fit-objekt použít funkci
# plot()
plot(lmfit2)

# The diagnostic plots show residuals in four different ways:
  
# Residuals vs Fitted. 
  # Used to check the linear relationship assumptions. 
  # A horizontal line, without distinct patterns is an indication 
  # for a linear relationship, what is good.
# Normal Q-Q. Used to examine whether the residuals are normally distributed. 
# Scale-Location (or Spread-Location). Used to check the homogeneity 
  # of variance of the residuals (homoscedasticity). Horizontal line with 
  # equally spread points is a good indication of homoscedasticity. 
# Residuals vs Leverage. Used to identify influential cases, that is extreme 
  # values that might influence the regression results when included or 
  # excluded from the analysis.

plot(lmfit2, 4, id.n = 10) # Cookovy vzdálenosti

# balíček car nabízí funkci residualPlots()
# Reziduální regresní grafy pro kontrolu linearity
  # Na ose X jsou jednotlivé prediktory
  # Na ose Y rezidua závislé proměnné po kontrole všech ostatních prediktorů
  # (tj. vyjma toho, který je na ose X)
car::residualPlots(lmfit2)

# funkce augment() z balíčků broom v případě lineární regrese doplní:
  # .fitted = predikované hodnoty
  # .resid = rezidua
  # .hat = leverage (páka)
  # .sigma = standardní odchylka reziduí po vyřazení daného případu
  # .cooksd = Cookovy vzdálenosti
  # .std.resid = standardizované reziduum
df_aug <- broom::augment(lmfit2, data = df_clean)
head(df_aug)

# Variance inflation factor
  # odhad toho, kolikanásobně jsou větší standardní chyby regresních
  # koeficientů oproti stavu, kdyby prediktory vůbec nekorelovaly
car::vif(lmfit2) %>% sqrt()

# Výběr pozorování s nejvyšší Cookovou vzdáleností
df_aug %>% 
  slice_max(order_by = .cooksd, n = 10)

# Výběr pozorování s nejvyšším reziduem
df_aug %>% 
  slice_max(order_by = abs(.resid), n = 10)

# Bodové grafy s vyznačenými vlivnými případy
df_aug %>% 
  ggplot(aes(vekr, deprese, color = .cooksd >= 0.01)) +
  geom_jitter(size = 3)

df_aug %>% 
  ggplot(aes(selfe, deprese, color = .cooksd >= 0.01)) +
  geom_jitter(size = 3)

# Srovnání dvou nezávislých skupin --------------------------

# * Nezávislý t-test ----------------------------------------
# Budeme testovat rozdíl v průměru známek mezi chlapci a dívkami

# Čištění dat
df_clean <- df %>% 
  drop_na(pohlavi, gpa)

### Grafy
# Graf hustoty pravděpodobnosti
df_clean %>% 
  ggplot(aes(gpa, color = pohlavi)) +
  geom_density(size = 2)

# Boxplot
df_clean %>% 
  ggplot(aes(pohlavi, gpa)) +
  geom_boxplot()

# Histogramy
df_clean %>% 
  ggplot(aes(gpa, y = stat(density))) +
  geom_histogram(binwidth = .25) +
  facet_wrap(~pohlavi)

# Q-Q grafy
df_clean %>% 
  ggplot(aes(sample = gpa)) +
  geom_qq() +
  geom_qq_line() +
  facet_wrap(~pohlavi)

# Deskriptivní statistiky
df_clean %>% 
  filter(!is.na(gpa)) %>% 
  group_by(pohlavi) %>% 
  summarise(
    n = n(),
    ci_mean(gpa, conf.level = .995),
    SD = sd(gpa),
    skew = psych::skew(gpa)
  )

# Shapiro-Wilk test of normality
# Není příliš užitečný, protože malý vzorek = malá statistická síla,
# velký vzorek = i drobné odchylky od normality budou signifikantní
df_clean %>% 
  group_by(pohlavi) %>% 
  rstatix::shapiro_test(gpa)

# Provedení t-testu
t_res <- t.test(gpa ~ pohlavi, data = df_clean, 
                conf.level = .995, # úroveň spolehlivosti
                var.equal = TRUE) # předpokládat shodu rozptylů?
t_res

# Velikost účinku, Cohenovo d
effectsize::cohens_d(gpa ~ pohlavi, data = df_clean)
effectsize::cohens_d(gpa ~ pohlavi, 
                     data = df_clean,
                     pooled_sd = FALSE,
                     ci = .995)

# Hednesovo g (v podstatě Cohenovo d s korekcí biasu)
  # u malých vzorků totiž má Cohenovo d sklon nadhodnocovat velikost účinku
effectsize::hedges_g(gpa ~ pohlavi, data = df_clean)


# * Mann–Whitney–Wilcoxon ---------------------------------
# Budeme testovat rozdíly mezi chlapci a dívkami ve známce z češtiny
df_clean <- df %>% 
  drop_na(pohlavi, zn_cj)

# Hlavní kvanltiy se neliší, protože nemáme moc unikátních hodnot 
df_clean %>% 
  group_by(pohlavi) %>% 
  summarise(
    quantiles(zn_cj)
)

# Ale z četností to vypadá, 
# že dívky dostávají častěji jedničky a dvojky než chlapci
# a chlpaci častěji trojky a čtyřky než dívky
df_clean %>% 
  filter(!is.na(zn_cj)) %>% 
  count(pohlavi, zn_cj) %>% 
  group_by(pohlavi) %>% 
  mutate(p = n/sum(n) *100)

# Provedení testu samotného
wilcox.test(zn_cj ~ pohlavi, data = df)


# Velikost účinku: Pravděpodobnost superiority
# musíme si připravit matici
m <- df_clean %>% 
  select(zn_cj, pohlavi) %>% 
  mutate(pohlavi = as.integer(pohlavi)) %>% 
  as.matrix()
head(m) # data musí mít formát matice, kde první sloupce je závislá proměnná
# a druhý sloupec nezávislá proměnná
# jsou povoleny jen numerické hodnoty
m %>% 
  RProbSup::A(
    design = 1, # mezisubjektový design je 1
    conf.level = .995 # úroveń spolehlivosti
    )
  
# Srovnání dvou závislých skupin ------------------------------------------

# * Párový t-test ---------------------------------------------------------
# Budeme testovat rozdíl mezi vnímanou vřelostí matky a otce

# Příprava dat
df_clean <- df %>% 
  drop_na(warm_m, warm_o)

# Převod na delší formát
df_long <- df_clean %>%
  select(id, warm_m, warm_o) %>% 
  mutate(warm_diff = warm_m - warm_o) %>% # vypočteme si i rozdíl
  pivot_longer(
    cols = c(warm_m, warm_o, warm_diff),
    names_to = "warm",
    values_to = "value"
  )
head(df_long)

# Deskriptivní statistiky
df_long %>% 
  group_by(warm) %>% 
  summarise(
    n = n(), 
    ci_mean(value, .995),
    SD = sd(value)
  )

# Zajímá nás obvykle i korelace mezi párem proměnných
cor.test(df_clean$warm_m, df_clean$warm_o, conf.level = .995)

# Grafy
# Histogram s křivkou normálního rozdělení
df_long %>% 
  filter(warm == "warm_diff") %>% 
  ggplot(aes(value)) +
  geom_histogram(aes(y = stat(density))) +
  stat_function(fun = dnorm, args = list(mean = 0.151, sd = 0.540),
                color = "red", size = 2)
# Q-Q graf
df_long %>% 
  filter(warm == "warm_diff") %>% 
  ggplot(aes(sample = value)) +
  stat_qq() +
  stat_qq_line()

# Provedení párového t-testu
t.test(df_clean$warm_m, df_clean$warm_o, 
       paired = TRUE, conf.level = .995)

# Výpočet Cohenova d
effectsize::cohens_d(df_clean$warm_m, df_clean$warm_o,
                     paired = TRUE,
                     ci = .995)



# * Wilcoxonův párový (sign-rank) test -------------------------------------
# Budeme opět testovat rozdíl mezi vnímanou vřelostí matky a otce

# Kvartily pro porovnání
df_long %>% 
  group_by(warm) %>% 
  summarise(
    quantiles(value)
  )

# Provedení testu
wilcox.test(df_clean$warm_m, df_clean$warm_o)

# Výpočet pravděpodobnosti superiority
df_clean %>% 
  select(warm_m, warm_o) %>% 
  as.matrix() %>% 
  RProbSup::A(
    design = 2, # vnitrosubjektový desing
    conf.level = .995 # úroveň spolehlivosti
  )

# Srovnání více nezávislých skupin ----------------------

# * One-way ANOVA ---------------------------------------
# Budeme testovat vztah mezi očekávaných vzděláním a průměrem známek

# Četnosti skupin
df %>% 
  count(ocek_vzd)

# Příprava dat
df_clean <- df %>% 
  select(id, skola, trida, ocek_vzd, gpa) %>% 
  na.omit()

# Boxplot
qplot(ocek_vzd, gpa, geom = "boxplot", data = df_clean)

# Průměry s intervaly spolehlivosti
df_clean %>% 
  ggplot(aes(ocek_vzd, gpa)) +
  stat_summary(fun.data = mean_cl_normal, 
               fun.args = list(conf = .995),
               geom = "errorbar", size = 1, width = 0.25) +
  stat_summary(fun = mean, geom = "point", size = 4)

# Histogram rozdělení gpa pro jednotlivé skupiny
df_clean %>% 
  ggplot(aes(gpa)) +
  geom_histogram(binwidth = .5) +
  facet_wrap(~ocek_vzd, scales = "free_y")

# Grafy hustoty pravděpodobnosti
df_clean %>% 
  ggplot(aes(gpa)) +
  geom_density() +
  facet_wrap(~ocek_vzd, scales = "free_y")

# Q-Q grafy
df_clean %>% 
  ggplot(aes(sample = gpa)) +
  geom_qq() +
  geom_qq_line() +
  facet_wrap(~ocek_vzd)

# Deskriptivní statistiky
df_clean %>% 
  group_by(ocek_vzd) %>% 
  summarise(
    n = n(),
    ci_mean(gpa),
    SD = sd(gpa),
    var = var(gpa),
    skew = skew(gpa),
    kurt = kurtosi(gpa)
  )

# Levenův test shody rozptylů
# Stejné problémy jako Shapiro-Wilk test
df_clean %>% 
  rstatix::levene_test(gpa ~ ocek_vzd, center = "mean")

# One-waz ANOVA
anova1 <- aov(gpa ~ ocek_vzd, data = df_clean)
summary(anova1)

# Welchova ANOVA
anova2 <- oneway.test(gpa ~ ocek_vzd, data = df_clean)
anova2

# Diagnostické grafy (stejné jako při použití lm())
plot(anova1)

# Velikost účinku omega2

effectsize::omega_squared(anova1,
                          alternative = "two.sided",
                          ci = .995)

### Kontrasty
levels(df_clean$ocek_vzd)

K <- matrix(
  c(-1/2,-1/2, 1/4, 1/4, 1/4, 1/4, # úrovně bez maturity vs s maturitou
       0,   0,-1/3,-1/3,-1/3, 1), # úrovně s maturitou (ale bez VŠ) vs VŠ 
  nrow = 2, byrow = TRUE,
  dimnames = list(
    c("bez/sMaturitou", "bez/sVŠ"),
    levels(df_clean$ocek_vzd)
  )

)
K

# Odhad samotných kontrastů s použitím balíčku multcomp
mc <-  multcomp::glht(anova1, linfct = K)
summary(mc)
confint(mc) # Intervaly spolehlivosti

# Post-hoc testy
# Tukeyho HSD má dobrou statistickou sílu, ale asi nejvíce předpokladů
# včetně normálního rozdělení a shody rozptylů v rámci skupin
TukeyHSD(anova1, conf.level = .995)
df_clean %>% 
  rstatix::tukey_hsd(gpa ~ ocek_vzd)

# games_howell test shodu rozptylů nepředpokládá,
# ale má menší statistickou sílu
df_clean %>% 
  rstatix::games_howell_test(gpa ~ ocek_vzd)

# * Kruskal-Wallisův Test ---------------------------------
# Budeme opět testovat vztah mezi očekávaných vzděláním a průměrem známek
kruskal.test(gpa ~ ocek_vzd, data = df_clean)

# Velikost účinku etaH2
# Jde vlastně také o odhad vysvětleného rozptylu (v pořadí/ranks)
df_clean %>% 
  rstatix::kruskal_effsize(gpa ~ ocek_vzd,
                           ci = FALSE, # vypnuto, ať se nezdržujeme
                           conf.level = .995,
                           nboot = 2000)


# Post-hoc srovnání pomocí Dunnova testu
# Neparametrická obdoba Tukeyho HSD testu
df_clean %>% 
  rstatix::dunn_test(gpa ~ ocek_vzd,
                     p.adjust.method = "holm")

# Vztah mezi dvěma kategorickými proměnnými -------------------------------

# Příprava dat
df_clean <- df %>% 
  drop_na(zn_cj, zn_mat)

# Custom funkce crosstabs() pro kontingenční tabulku
crosstabs(df, "zn_cj", "zn_mat")

# Kontingenční tabulky
# Funkce CrossTable() z balíčku gmodels
gmodels::CrossTable(
  df$zn_cj, df$zn_aj,
  prop.r = TRUE, prop.c = FALSE, prop.t = FALSE,
  prop.chisq = FALSE,
  expected = TRUE, 
  asresid = TRUE,
  format = "SPSS")

# Base funkce chisq.test()
chisq.test(df$zn_cj, df$zn_mat)
chisq.test(df$zn_cj, df$zn_mat,
           simulate.p.value = TRUE, # p-hodnoty na základě simulace
           B = 1e4) # počet simulací

# Velikost pčinku, Cramérovo V
effectsize::cramers_v(df$zn_cj, df$zn_mat,
                      ci = .995, # hladina spolehlivosti
                      alternative = "two.sided", # typ intervalu
                      adjust = TRUE) # korekce bias

# Mosaic plot
df_clean %>% 
  plot_mosaic(xvar = "zn_cj", yvar = "zn_mat") +
  scale_fill_manual(values = palette.colors(n = 5, palette = "R4"))

# Using geom_count()
df_clean %>% 
  ggplot(aes(zn_cj, zn_mat)) + 
  geom_count() +
  scale_size_continuous(
    range = c(1, 10)
  )




