# Instalace balíčků, pokud nejsou instalovány
pckg <- c("psych", "broom", "car", "effectsize",
          "corrplot", "rstatix", "lm.beta",
          "tidyverse", "ggmosaic", "GGally",
          "apaTables", "gmodels")

pckg <- pckg[!(pckg %in% installed.packages()[,"Package"])]
if(length(pckg)) install.packages(pckg)

library(tidyverse)
library(ggmosaic)
library(GGally)

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

# Přehled dat
glimpse(df)

ci_mean <- function(x, level = .95, na.rm = TRUE) {
  if (na.rm == TRUE) {
    x <- x[!is.na(x)]
  }
  
  n <- length(x)
  df <- n - 1
  s <- sd(x)
  m <- mean(x)
  
  se <- s/sqrt(n)
  alpha <- 1 - .95
  t <- qt(1 - alpha/2, df = df)
  me <- t*se
  ci_lower <- m - me
  ci_upper <- m + me
  
  tibble(m = m,
         ci_lower = ci_lower,
         ci_upper = ci_upper,
         sd = s,
         se = se,
         n = n)
  
}

quartiles <- function(x,
                      na.rm = TRUE) {
  tibble(
    stat = c("Min", "q25", "Mdn", "q75", "Max"),
    value = quantile(x, 
                     probs = c(0, .25, .5, .75, 1), 
                     na.rm = na.rm)
  ) %>% 
    pivot_wider(names_from = stat, values_from = value)
}


# Korelace --------------------------------------------
# Co korelace, to scatterplot
df %>% 
  ggplot(aes(optim, ziv_sp)) +
  geom_jitter(alpha = 0.25, size = 3,
              width = 0.5, height = 0.1) +
  geom_smooth(method = "lm", color = "blue", 
              se = FALSE, linewidth = 3) +
  geom_smooth(color = "red", fill = "red",
              alpha = .1)

# Výběr proměnných, mezi kterými chceme zjistit korelace
variables <- c("optim", "ziv_sp", "selfe", "effi", "zdravi", "deprese")
df_cor <- df %>% select(all_of(variables))
df_cor

# Matice grafů
df_cor %>% 
  GGally::ggpairs()

df_cor %>% 
  GGally::ggpairs(
    upper = list(continuous = wrap("cor")),
    diag = list(continuous = wrap("barDiag")), 
    lower = list(continuous = wrap("smooth_loess", alpha = 0.1,
                                   position = "jitter"))
  )
# https://ggobi.github.io/ggally/articles/ggally_plots.html


# cor() funkce base R
df_cor %>% 
  cor(use = "pairwise") %>% 
  round(2)

df_cor %>% 
  cor(use = "pairwise",
      method = "spear") %>% 
  round(2)

df_cor %>% 
  cor(use = "pairwise",
      method = "kendall") %>% 
  round(2)

# Balíček psych
r <- psych::corr.test(df_cor,
                      use = "pairwise", # Způsob zacházení s chybějícími daty
                      adjust = "holm", # jakou korekci p-hodnot použít
                      method = "spearman", # typ korelace
                      alpha = 0.005) # hladina spolehlivosti bude 1 - alpha

print(r, short = FALSE)
r$stars
str(r)

# Bootstrap
psych::cor.ci(df_cor, 
              n.iter = 2000)

# Balíček corrplot
r <- cor(df_cor, use = "pairwise")

corrplot::corrplot.mixed(r)
corrplot::corrplot.mixed(r, 
                         lower.col = "black",
                         order = "hclust") # Seřadit na základě klustrové analýzy
df_cor %>% 
  apaTables::apa.cor.table()

# Lineární regrese ----------------------------------------------
# V rámci cvičné analýzy budeme predikovat depresivitu na základě 
# pohlaví, věku, průměru známek, důvěry v rodiče, důvěry ve vrstevníky, 
# sebehodnocení a optimismu.

# Vyčištění dat
df_clean <- df %>% 
  select(id, skola, trida, deprese, pohlavi, vekr, 
         gpa, duv_v, duv_r, selfe, optim) %>% 
  drop_na()


df_clean %>% 
  select(where(is.numeric), -c(id:trida)) %>% 
  cor() %>% 
  corrplot::corrplot.mixed(lower = "number", # Hodnoty korelací dole
                           upper = "circle", # Nahoře elipsy
                           lower.col = "black") # Seřadit na základě klustrové analýzy

df_clean %>% 
  ggplot(aes(pohlavi, deprese)) +
  geom_boxplot()

# Všechny bivariační vztahy
df_clean %>% 
  select(pohlavi:optim, deprese) %>% 
  GGally::ggpairs(
    upper = list(continuous = wrap("cor"),
                 combo = wrap("facethist")),
    diag = list(continuous = wrap("barDiag")),
    lower = list(continuous = wrap("smooth_loess", alpha = 0.05,
                                   position = "jitter"),
                 combo = wrap("box"))
  )

# Samotný odhad modelu
lmfit1 <- lm(deprese ~ pohlavi + vekr, 
             data = df_clean) 
lmfit1
summary(lmfit1)

lmfit2 <- lm(deprese ~ pohlavi + vekr + 
               gpa + duv_v + duv_r + selfe + optim, 
             data = df_clean)
summary(lmfit2)

lmfit2 <- update(lmfit1, # Původní model
                 . ~ . + gpa + duv_v + duv_r + selfe + optim)
summary(lmfit2)

lmfit2 %>% 
  lm.beta::lm.beta() %>% 
  summary()

apaTables::apa.reg.table(lmfit2)
apaTables::apa.reg.boot.table(lmfit2)
apaTables::apa.reg.table(lmfit1, lmfit2)

# Funkce pro výpočet rozdílu v R2 mezi modely
r2_diff <- function(fit1, fit2) {
  c(R2_diff = summary(fit2)$r.squared - summary(fit1)$r.squared,
    R2_adj_diff = summary(fit2)$adj.r.squared - summary(fit1)$adj.r.squared)
}

r2_diff(lmfit1, lmfit2)
anova(lmfit1, lmfit2)

# Koeficienty modelu v přehledné tabulce (tibble)
broom::tidy(lmfit2, conf.int = TRUE, conf.leve = 0.99)

# Informace o fitu modelu
broom::glance(lmfit1)
broom::glance(lmfit2)

bind_rows(lmfit1 = broom::glance(lmfit1), 
          lmfit2 = broom::glance(lmfit2),
          .id = "fit")

# Informací o jednotlivých případech
df_aug <-  broom::augment(lmfit2, se_fit = TRUE, 
                          interval = "confidence")

df_aug <- df_aug %>% 
  mutate(id = row_number()) %>% 
  select(id, everything())

df_aug

# .fitted = predikované hodnoty
# .lower a .upper = interval spolehlivosti nebo predikční interval (jeho meze)
# .se.fit	- stanradní chyba predikovaných hodnot
# .resid	= reziduum
# .hat	= páka
# .sigma = stanardní chyba reziduí po vyřazení daného případu
# .std.resid = standardizované reziduum

# Diagnostika

plot(lmfit2)

plot(lmfit2, 4, 
     id.n = 10) # vyznačit 10 případů s největší Cookovou vzdáleností

# Rezidua vs hodnoty prediktorů
car::residualPlots(lmfit2)
car::ceresPlots(lmfit2)

# Added variable plots
car::avPlots(lmfit2)

# Vlivné či extrémní případy
car::influencePlot(lmfit2)
car::influenceIndexPlot(lmfit2)

# qq-graf reziduí
car::qqPlot(lmfit2, main="QQ Plot")

# Variance Infation Factor
car::vif(lmfit2) 

# Durbin-Watson test
car::durbinWatsonTest(lmfit2)


# S "augmentovanými" daty můžeme dělat cokoli
df_aug %>% 
  slice_max(order_by = .cooksd, n = 10)

df_aug %>% 
  slice_max(order_by = abs(.resid), n = 10)

df_aug %>% 
  ggplot(aes(.hat, .std.resid, size = .cooksd)) +
  geom_text(aes(label = id))


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

# Nezávislý t-test -------------------------------------------
df_clean <- df %>% 
  drop_na(pohlavi, gpa)

# Histogram s relativní četností
# Pozn. hustato pravděpodobnosti * binwidth = relativní četnost
df_clean %>% 
  ggplot(aes(gpa, y = after_stat(density)*0.25)) + # Hustota pravděpodobnosti 
  geom_histogram(binwidth = 0.25) +
  facet_wrap(~pohlavi)

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

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

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

df_clean %>% 
  group_by(pohlavi) %>% 
  summarise(
    n = n(), # počet případů
    ci_mean(gpa, level = .95), # průměr s intervalem spolehlivosti
    skew = psych::skew(gpa) # šikmost
  )

df_clean %>% 
  group_by(pohlavi) %>% 
  rstatix::shapiro_test(gpa)

df_clean %>% 
  rstatix::levene_test(gpa ~ pohlavi)

t.test(gpa ~ pohlavi, # Závislá proměnná ~ nezávislá proměnná
       data = df_clean, # Použitý dataset
       conf.level = .995, # úroveň spolehlivosti
       var.equal = TRUE) # předpokládat shodu rozptylů?

df_clean %>% 
  rstatix::t_test(gpa ~ pohlavi,
                  conf.level = 0.995,
                  var.equal = TRUE,
                  detailed = TRUE)


df_clean %>% 
  ggplot(aes(pohlavi, gpa)) +
  geom_pointrange(stat = "summary",
                  fun.data = "mean_cl_normal") +
  geom_line(stat = "summary",
            aes(group = 1))

effectsize::cohens_d(gpa ~ pohlavi, 
                     data = df_clean,
                     ci = 0.995)


# * Wilcoxon-Mann–Whitney (WMW) U-test --------------------------------------
df_clean <- df %>% 
  drop_na(pohlavi, zn_cj)

df_clean %>% 
  group_by(pohlavi) %>% 
  summarise(
    quartiles(zn_cj)
  )

freq <-  df_clean %>% 
  count(pohlavi, zn_cj) %>% 
  group_by(pohlavi) %>% 
  mutate(p = n/sum(n))
freq


freq %>% 
  ggplot(aes(zn_cj, p, fill = pohlavi)) +
  geom_col(position = "dodge")

freq %>% 
  mutate(zn_cj = factor(zn_cj)) %>% 
  ggplot(aes(pohlavi, p, fill = fct_rev(zn_cj))) +
  geom_col()


wilcox.test(zn_cj ~ pohlavi, data = df_clean)

df_clean %>% 
  rstatix::wilcox_test(zn_cj ~ pohlavi)

effectsize::p_superiority(zn_cj ~ pohlavi, 
                          data = df_clean,
                          iterations = 2000, # Bootstrap s 2000 iteracemi
                          parametric = FALSE) # Nepředpokládat normální rozdělení

# Párové testy -------------------------------------------------

# * Explorace dat ---------------------------------------------
df_clean <- df %>% 
  drop_na(warm_m, warm_o) %>% 
  mutate(warm_diff = warm_m - warm_o)

df_clean %>% 
  summarise(
    ci_mean(warm_diff),
    quartiles(warm_diff)
  )

df_long <- df_clean%>% 
  pivot_longer(starts_with("warm"),
               names_to = "parent",
               values_to = "warm") 

df_long %>% 
  group_by(parent) %>% 
  summarise(ci_mean(warm),
            quartiles(warm),
            skew = psych::skew(warm),
            kurt = psych::kurtosi(warm))
  
df_clean %>% 
  select(starts_with("warm")) %>% 
  psych::cor.ci(n.iter = 2000)

df_clean %>% 
  ggplot(aes(warm_m, warm_o)) +
  geom_jitter(width = 0.05, height = 0.05,
              alpha = 0.2) +
  coord_fixed() +
  geom_abline(slope = 1) +
  scale_x_continuous(limits = c(0.95, 4.05)) +
  scale_y_continuous(limits = c(0.95, 4.05))

# Histogram s křivkou normálního rozdělení
df_clean %>% 
  ggplot(aes(warm_diff)) +
  geom_histogram(aes(y = stat(density))) +
  stat_function(fun = dnorm, 
                args = list(mean = mean(df_clean$warm_diff), 
                            sd = sd(df_clean$warm_diff)),
                color = "red", size = 2)

# Q-Q graf
df_clean %>% 
  ggplot(aes(sample = warm_diff)) +
  stat_qq() +
  stat_qq_line()

# * Párový t-test --------------------------------------------
out <- t.test(df_clean$warm_m, df_clean$warm_o, # Porovnávaný pár proměnných
              paired = TRUE, conf.level = .995)
broom::tidy(out)

# Cohenovo d
effectsize::cohens_d(df_clean$warm_m, df_clean$warm_o,
                     paired = TRUE,
                     ci = .995)

# * Wilcoxonův párový (sign-rank) test --------------------------------
wilcox.test(df_clean$warm_m, df_clean$warm_o,
            paired = TRUE)

# Velikost účinku
effectsize::p_superiority(df_clean$warm_m, df_clean$warm_o, 
                          data = df_clean,
                          paired = TRUE,
                          iterations = 2000, # Bootstrap s 2000 iteracemi
                          parametric = FALSE) # Nepředpokládat normální rozdělení

# Vztah mezi dvěma kategorickými proměnnými ------------------------------
df_clean <- df %>% 
  drop_na(zn_cj, zn_mat) %>% 
  mutate(across(starts_with("zn_"), factor))

gmodels::CrossTable(df_clean$zn_cj, df_clean$zn_mat,
                    prop.chisq = FALSE,
                    asresid = TRUE,
                    digits = 2,
                    format = "SPSS",
                    chisq = TRUE)

Xsq <- chisq.test(df_clean$zn_cj, df_clean$zn_mat)
Xsq <- chisq.test(df_clean$zn_cj, df_clean$zn_mat,
                  simulate.p.value = TRUE,
                  B = 2000)
Xsq
Xsq$observed
Xsq$expected
Xsq$residuals
Xsq$stdres

df_clean %>% 
  ggplot() +
  ggmosaic::geom_mosaic(aes(x = product(zn_mat, zn_cj), fill = zn_mat))

df_clean %>% 
  ggplot() +
  ggmosaic::geom_mosaic(aes(x = product(zn_cj, zn_mat), fill = zn_cj))