# Načtení balíčků a import dat --------------------------------------------

# Balíčky
library(tidyverse)
library(skimr)
source("https://is.muni.cz/el/fss/podzim2022/PSYn5320/um/r101_funs.R")

# Import dat
df <- readRDS(
  url("https://is.muni.cz/el/fss/podzim2022/PSYn5320/um/datasets/eukids.rds", 
      "rb"))

# Přehled dat
glimpse(df)

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

# Custom skim funkce
my_skim <- skim_with(
  numeric = sfl(
    Skw = ~psych::skew(.),
    Krt = ~psych::kurtosi(.)
  )
)

# Lineární regrese --------------------------------------------------------
# Predikujte životní spokojenost (ziv_sp) na základě: věku 
# (vekr), pohlaví (pohlavi), optimismu (optim), monitorování (monit), vztahu s 
# rodiči (vztahsr) a důvěrnosti s vrstevníky (duv_v) a průměru známek (gpa). 
# Modely odhadněte ve dvou krocích. V prvním kroku vložte „nepsychologické“ 
# proměnné (věk, pohlaví a průměr známek), v druhém kroku všechny ostatní 
# prediktory. Ověřte statistické předpoklady regrese (linearitu, normální 
# rozdělení reziduí, nezávislost reziduí, absenci extrémních/vlivných případů, 
# absenci vysoké multikolinearity). Srovnejte oba modely z hlediska shody s daty 
# (vysvětleného rozptylu). Interpretujte lépe fitující model: celkový vysvětlený 
# rozptyl a efekty jednotlivých prediktorů.

df_clean <- df %>% 
  select(id, ziv_sp, pohlavi, vekr, gpa, monit, vztahsr, duv_v, optim) %>% 
  drop_na()

fit1 <- lm(ziv_sp ~ pohlavi + vekr + gpa, data = df_clean)
summary(fit1)
lm.beta::lm.beta(fit1)$standardized.coefficients

fit2 <- lm(ziv_sp ~ pohlavi + vekr + gpa +
             monit + vztahsr + duv_v + optim, data = df_clean)
summary(fit2)
lm.beta::lm.beta(fit2)$standardized.coefficients

r2_diff(fit1, fit2)
car::avPlots(fit2)
car::residualPlots(fit2)
plot(fit2)

broom::augment(fit2) %>% 
  ggplot(aes(.std.resid)) +
  geom_histogram()
  
broom::augment(fit2) %>% 
  ggplot(aes(sample = .std.resid)) +
  geom_qq() +
  geom_qq_line()

car::durbinWatsonTest(fit2)
car::vif(fit2)
r <- df_clean %>% 
  mutate(pohlavi = as.integer(pohlavi)) %>% 
  select(-id) %>% 
  cor()

r %>% 
  corrplot::corrplot.mixed(lower = "number", # Hodnoty korelací dole
                           upper = "ellipse", # Nahoře elipsy
                           lower.col = "black") # Seřadit na základě klustrové analýzy
  
fit3 <- lm(ziv_sp ~ pohlavi + vekr + gpa +
             vztahsr + duv_v + optim, data = df_clean)
summary(fit3)
car::vif(fit3)
lm.beta::lm.beta(fit3)$standardized.coefficients

df_clean %>% 
  mutate(pohlavi = as.integer(pohlavi)) %>% 
  select(-id) %>% 
  ppcor::spcor() %>% 
  .$estimate %>% 
  t() %>% 
  round(digits = 2)

# Korelace ----------------------------------------------------------------
# Pomocí korelační analýzy prozkoumejte, které negativní pocity (neg1, neg2, 
# neg3) nejtěsněji souvisejí s depresivitou (deprese). Vytvořte bodové graf 
# těchto vztahů.

df_clean <- df %>% 
  select(id, deprese, neg1:neg3) %>% 
  drop_na()

df_clean <- df_clean %>% 
  mutate(across(neg1:neg3, ~.x * -1))

# neg1 - Negativní pocity - osamělost, strach, vina
# neg2 - Negativní pocity - unava, stres
# neg3 - Negativní pocity - agresivní

plot_scatter_mult(df_clean, 
                  xvar = c("neg1", "neg2", "neg3"),
                  yvar = "deprese",
                  geom = "jitter")
df_clean %>% 
  select(-id) %>% 
  psych::corr.test() %>% 
  print(short = FALSE)

psych::r.test(n = 765,
              r12 = 0.57, 
              r13 = 0.41, 
              r23 = 0.44)

# Srovnání dvou skupin ---------------------------------------------------
# Vyberte si jen jednu kohortu (kohorta) a pomocí nezávislého t-testu srovnejte 
# chlapce a dívky v monitorování rodiči (monit). Ověřte předpoklad normálního 
# rozdělení v rámci skupin a předpoklad shody rozptylů. Vypočtěte relevantní 
# deskriptivní statistiky pro obě skupiny. Vypočtěte velikost účinku společně 
# s intervalem spolehlivosti. Znázorněte rozdíly mezi skupinami graficky 
# (např. boxploty nebo grafy průměrů s chybovými úsečkami pro intervaly 
# spolehlivosti)

df_clean <- df %>% 
  select(kohorta, pohlavi, monit) %>% 
  drop_na() %>% 
  filter(kohorta == "starší kohorta")

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

df_clean %>% 
  ggplot(aes(monit, color = pohlavi)) +
  geom_density(size = 2)

df_clean %>% 
  ggplot(aes(monit, y = after_stat(density))) +
  geom_histogram(binwidth = 0.2) +
  facet_wrap(~pohlavi)

df_clean %>% 
  ggplot(aes(pohlavi, monit)) +
  geom_pointrange(stat = "summary",
                  fun.data = mean_cl_normal)

df_clean %>% 
  group_by(pohlavi) %>% 
  summarise(
    ci_mean(monit),
    SD = sd(monit),
    n = n()
  )

t.test(monit ~ pohlavi, data = df_clean, var.equal = TRUE)

effectsize::cohens_d(monit ~ pohlavi,
                     data = df_clean)

# Dvě závislá měření ------------------------------------------------------
# Vyberte si jen jedno pohlaví (pohlavi) a kohortu (kohorta), např. pouze 
# mladší dívky, a pomocí párového t-testu srovnejte důvěrnost s rodiči a 
# vrstevníky (duv_r a duv_v). Ověřte předpoklad normálního rozdělení rozdílových 
# skórů. Vypočtěte relevantní deskriptivní statistiky. Vypočtěte velikost 
# účinku společně s intervalem spolehlivosti. Znázorněte rozdíly mezi důvěrností 
# s rodiči a vrstevníky graficky (např. boxploty nebo grafy průměrů s chybovými 
# úsečkami pro intervaly spolehlivosti).

df_clean <- df %>% 
  filter(pohlavi == "zenske", kohorta == "mladší kohorta") %>% 
  select(id, pohlavi, kohorta, duv_r, duv_v) %>% 
  drop_na()

df_clean <- df_clean %>% 
  mutate(duv_diff = duv_r - duv_v, # rozdíl mezi vřelostí matky a otce
         duv_tot = (duv_r + duv_v)/2, # průměrná vřelost žáka pro oba rodiče
         grand_mean = mean(duv_tot), # celkový průměr obou typů vřelosti
         btw_eff = duv_tot - grand_mean, # odhad efektu žáka
         adj_duv_r = duv_r - btw_eff,
         adj_duv_v = duv_v - btw_eff
  )

df_clean %>% 
  ggplot(aes(duv_diff)) +
  geom_histogram()

df_clean %>% 
  my_skim(duv_r, duv_v, duv_diff)

df_clean %>% 
  ggplot(aes(sample = duv_diff)) +
  geom_qq()+
  geom_qq_line()

df_clean %>% 
  ggplot(aes(duv_r, duv_v)) +
  geom_jitter(alpha = 0.25) +
  geom_smooth()

df_clean %>% 
  select(starts_with("adj_")) %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(name, value)) +
  geom_pointrange(stat = "summary",
                  fun.data = mean_cl_normal)

t.test(df_clean$duv_r, df_clean$duv_v, 
       paired = TRUE)

effectsize::cohens_d(df_clean$duv_r, 
                     df_clean$duv_v,
                     paired = TRUE)