# Balíčky

# Instalovat ještě nenainstalované balíčky
list.of.packages <- c("mice", "tidyverse", "naniar", "haven",
                      "simputation", "VIM")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

# Načíst balíčky
library(mice)
library(tidyverse)
library(naniar)
library(simputation)
library(VIM)



# Základní funkce -----------------------------------------
x <- c(NA, NaN, Inf, ".", "missing")
x

# Které hodnoty jsou NA?
is.na(x)

# Je aspoň jedna hodnota NA?
naniar::any_na(x)

# Jsou všechny hodnoty NA?
naniar::all_na(x)
naniar::all_miss(x)

# Jsou všechny hodnoty kompletní?
naniar::all_complete(x)

# Import dat -----------------------------------------------------------------
# Použitý dataset pochází z výkumu, který zkoumal souběžnou validitu 
# nového testu kognitivních schopností (LEM) u romských a neromských dětí
df <- haven::read_sav("data/lem.sav") %>% 
  mutate(
    across(c(sch_type, school, class, sex, ethnic, 
             starts_with("edu_")), as_factor),
    b_total = b1 + b2
  )

glimpse(df)
summary(df)

# Celkový počet NA
naniar::n_miss(df)

# Celkový podíl NA
naniar::prop_miss(df)

# Celkový podíl validních hodnot
naniar::prop_complete(df)

# Explorace chybějících dat ------------------------------

# ** Tabulky ------------------------------------------------
# Počet a podíl chybějících hodnot v jednotlivých proměnných
miss_var_summary(df)

# Otázka: Ve kterých proměnných chybí nejvíce hodnot?
  
  
# Kolik hodnot chybí kolika proměnným
miss_var_table(df)

# Srovnání etnik
df %>% 
  group_by(ethnic) %>% # Rozdělit podle ethnic
  miss_var_summary() %>% # Vypočíst počet a podíl chybějících hodnot
  pivot_wider(names_from = ethnic, values_from = c(n_miss, pct_miss)) %>% # Převést na wider formát
  mutate(diff = pct_miss_nerom  - pct_miss_Rom) %>%  # Vypočíst rozdíl v podílu NA
  arrange(-abs(diff)) # seřadit podle vypočteného rozdílu

# Otázka: Ve kterých se obě skupiny nejvíce liší, 
# co se týče podílu chybějících hodnot? Jaké vás napadá vysvětlení?
  
  

# počet chybějících hodnot jednotlivých respondentů
miss_case_summary(df)

miss_case_summary(df) %>% 
  filter(pct_miss > .5)

# kolika respondentům chybí kolik hodnot
miss_case_table(df)

miss_case_table(df) %>% 
  ggplot(aes(n_miss_in_case, pct_cases)) +
  geom_col()

# Případy s největším počtem NA
df %>% 
  slice(12:17) %>% 
  print(width = Inf)

# Otázka: V čem se shodují žáci s nejvtším počtem chybějících hodnot

# Funkcí add_n_miss() můžete přidat do datasetu sloupec 
# s počtem chybějících hodnot v daném řádku. 

df <- df %>% 
  add_n_miss()

df %>% 
  select(n_miss_all, everything()) %>% 
  arrange(desc(n_miss_all))

df %>% 
  add_n_miss(starts_with("lem"), label = "lem_miss") %>% 
  relocate(lem_miss_vars)

# ** Grafy ------------------------------------------------------
# zobrazí počet nebo podíl chybějících hodnot v jendotlivých proměnných.
gg_miss_var(df)

# Podíl chybějících hodnot podle skupin
gg_miss_var(df, facet = ethnic, show_pct = TRUE)

gg_miss_var(df, facet = sch_type, show_pct = TRUE)

gg_miss_var(df, facet = sex, show_pct = TRUE)

# V jakých proměnných se z hlediska podílu chybějících hodnot nejvíce liší 
# Romové vs neromové a žáci ze zvláštní vs základní školy?

# Které proměnné mají často společně missing values
vis_miss(df, cluster = TRUE) + 
  coord_flip() # prohodit osy

# Nejčastější patterns of missing values
gg_miss_upset(df)

gg_miss_upset(df, 
              nsets = 10, # maximální počet proměnných
              nintersects = 10) # maximální počet vzorců

# Hodnoty kterých proměnných mají tendenci chybět společně?

# Srovnat více skupin v podílu chybějících hodnot umožňuje také heatmapa.
gg_miss_fct(df, fct = ethnic)

gg_miss_fct(df, fct = school)

gg_miss_fct(df, fct = sch_type)

# Která škola school (kvůli anonymizaci známe jen její numerické id) 
# se od ostatních nejvíce odlišuje z hlediska podílu chybějících hodnot 
# a o které proměnné se jedná?

# Scatterploty s chybějícími hodnotami můžeme získat 
# pomocí funkce geom_miss_point()
ggplot(df, aes(raven, lem_total)) + 
  geom_miss_point(position = position_jitter(width = 0.5))

ggplot(df, aes(maths, lem_total)) + 
  geom_miss_point(position = position_jitter(width = 0.5))

# Vypadá to, že by chybění v lem_total záviselo na hodnotách raven nebo 
# maths nebo naopak chybění v raven nebo maths záviselo na hodnotě lem_total?

# Nabular data ----------------------------------------------------
# Tzv. shadow matrix obsahuje stejný počet řádků a sloupců jako původní 
# dataset, ale místo původních hodnot obsahují pouze informace o tom, 
# zda je nějaká hodnota chybějící (NA), nebo není (!NA). 
# Nabular data pak vzniknout připojením shadow matrix k původnímu 
# datasetu jako nové sloupce.

# Pouze shadow matrix
as_shadow(df) 

# argumentem (only_miss = TRUE) vyřadíme ze shadow matrix proměnné bez chybějících hodnot
df_nabular <- bind_shadow(df, only_miss = TRUE)

glimpse(df_nabular)

# Shadow matrix pak můžeme využít pro porovnání skupin vymezených 
# podle chybějících dat v nějaké proměnné, abychom zjistili, 
# zda chybění nějak systematicky nesouvisí s jinou proměnnou, 
# jejíž hodnoty známe.

df_nabular %>% 
  group_by(raven_NA) %>% 
  summarise(n = sum(!is.na(lem_total)),
            M_lem = mean(lem_total, na.rm = TRUE),
            SD_lem = sd(lem_total, na.rm = TRUE))

df_nabular %>% 
  group_by(maths_NA) %>% 
  summarise(n = sum(!is.na(lem_total)),
            M_lem = mean(lem_total, na.rm = TRUE),
            SD_lem = sd(lem_total, na.rm = TRUE))

df_nabular %>% 
  group_by(edu_f_NA) %>% 
  summarise(n = sum(!is.na(lem_total)),
            M_lem = mean(lem_total, na.rm = TRUE),
            SD_lem = sd(lem_total, na.rm = TRUE))

# Otázka: Na chybějí ve které proměnné závisí průměrný skór lem_total 
# nejsilněji? Jak to může souviset s následujícícm vztahem?
df_nabular %>% 
  count(sch_type, edu_f_NA) %>% 
  group_by(sch_type) %>% 
  mutate(p = n/sum(n))

# Samozřejmě shadow matrix můžeme využít i k tvorbě grafů. 
# Zde je např. graf hustoty pravděpodobnosti lem_total podle 
# chybění v reading (barevně) a podle ethnic (fazety)
df_nabular %>% 
  ggplot(aes(lem_total, color = reading_NA)) +
  geom_density() +
  facet_wrap(~ethnic)

# Liší se nějak výrazně rozdělení lem_total v závislosti 
# na chybějí v proměnné reading u Romů nebo neromů?

# Imputace dat -----------------------------------------------------
# V analýze se zaměříme primárně na to, jaký efekt má škola (sch_type) 
# na zvládnutí základních prvků gramotnosti - zaměříme se na čtení (reading).

# Čtení bylo testováno školním testem. Předpokládáme, že v základní škole 
# se dítě naučí více než ve zvláštní. 
# V duchu projektu budeme zvažovat, že efekt typu školy 
# (tj. jak dobře dokáže škola základy žáka naučit) se liší podle etnicity (ethnic).

# I když bychom si to úplně nepřáli, je možné, že roli bude hrát také 
# pohlaví žáka (sex). 
# Konečně musíme zohlednit také intelekt (lem), 
# protože schopnost čtení je jím samozřejmě ovlivněna.

# ** Imputace průměru ----------------------
# Nejjednoduší metodou imputace je jako chybějící hodnoty imputovat 
# průměr dané proměnné (vypočtený na základě nechybějících hodnot). 
# Pro kategorické proměnné je imputace modu obdobou imputace průměru 
# u numerických proměnných.

df_imp_mean <- df %>% 
  bind_shadow() %>% 
  add_label_shadow() %>% 
  mutate(
    across(where(is.numeric), impute_mean)
  )

# Zásadní problém je ten, že imputované hodnoty jsou vlastně konstantou, 
# což samozřejmě může zásadně ovlivnit univariančí i multivariační rozdělení 
# našich dat (viz níže)
ggplot(df_imp_mean, aes(lem_total, reading, color = any_missing)) +
  geom_point()

ggplot(df_imp_mean, aes(reading, fill = reading_NA)) +
  geom_histogram()

# Kvůli nižší variabilitě se sníží vztahy s ostatními proměnnými.

# ** Hot-deck imputace -----------------------------------------

# Principem je jednoduše nahradit chybějící hodnoty poslední 
# nechybějící hodnotou. Pro doplnění “validnějších” hodnot má smysl 
# stanovit si tzv. “doménovou” proměnnou nebo proměnné, 
# aby se chybějící hodnoty doplňovaly pouze v rámci téže “domény”.

# Vypočtení "globálního" skóru inteligence
df <- df %>% 
  mutate(global = pick(raven, b1, b2, lem1:lem6) %>% 
           mutate(across(everything(), .fns = scale)) %>% 
           rowMeans(na.rm = TRUE))

# Samotná imputace
df_imp_hd <- df %>% 
  VIM::hotdeck(variable = c("lem_total", "reading"), # do kterých proměnných imputovat data
               domain_var = c("sch_type", "ethnic"), # vymezení "domén"
               ord_var = "global") # seřadit podle globálního skóru 

# Pak se opět můžeme podívat na vztah obou proměnných po imputaci. 
# Můžeme nyní zkusit použít funkci marginplot() z balíčku VIM.
# Draw a margin plot
df_imp_hd %>% 
  select(lem_total, reading, ends_with("_imp")) %>% 
  VIM::marginplot(delimiter = "_imp")

# Otázka: Jak si vede hot-deck imputace ve srovnání s mean imputation

# ** kNN imputace ----------------------------------------
# Doplnění hodnot na základě "nejbližších" (nejpodobnějších) případů
# (nearest neighbours)
# 

df_imp_knn <- df %>% 
  VIM::kNN(variable = c("lem_total", "reading"), # Imputované proměnné
           dist_var = c("sch_type", "ethnic", # proměnné pro výpočet vzdáleností
                        "raven", "b_total"),
           k = 5, # počet neighbors
           numFun = weighted.mean, # funkce pro výpočet hodnoty
           weightDist = TRUE) # Vážit vzdáleností

# Tohle vypadá jako artefakt
df_imp_knn %>% 
  select(lem_total, reading, ends_with("_imp")) %>% 
  marginplot(delimiter = "_imp")

random <- function(x) {
  sample(x, size = 1)
}

df_imp_knn <- df %>% 
  VIM::kNN(variable = c("lem_total", "reading"), # Imputované proměnné
           dist_var = c("sch_type", "ethnic", # proměnné pro výpočet vzdáleností
                        "raven", "b_total"),
           k = 5, # počet neighbors
           numFun = random, # funkce pro výpočet hodnoty
           weightDist = FALSE) # Vážit vzdáleností

df_imp_knn %>% 
  select(lem_total, reading, ends_with("_imp")) %>% 
  marginplot(delimiter = "_imp")

# ** Vícenásobná imputace (multiple imputation) -----------------
# Slabinou všech předchozích metod, které jsme si ukázali, je to, 
# že zanedbávají nejistotu imputovaných hodnot, ale pracují s nimi, 
# jako bych se jednalo o skutečně pozorované hodnoty 
# a my měli k dispozici kompletní data. 
# To obvykle vede k podhodnocení standardních chyb, 
# protože ty počítají pouze s výběrovou variabilitou, 
# nikoli imputační variabilitou – a tu je třeba zohlednit.

# Doporučení pro specifikaci imputačních modelů jsou následující:
# 1. Zahrňme všechny proměnné které budou součastní finálního modelu
#    (cílové proměnné)
# 2. Zahrňme všechny proměnné, které souvisejí s absencí odpovědí 
#    v cílových proměnných
# 3. Zahrňme proměnné vysvětlující významné podíl rozptylu v proměnných, 
#    které plánujem imputovat. 
# 4. Vyřaďte ty proměnné z kroku 2 a 3, 
#    které mají příliš chybějících hodnot
#    (hlavně v podskupině nekomplentích případů z cílových proměnných)

# Dejme tomu, že chceme imputovat hodnoty pěti numerických proměnných 
# a imputační modely specifikujeme takto:
form_list <- c(
  raven = raven ~ sex + sch_type*ethnic +  b_total + reading + maths + lem_total,
  b_total = b_total ~ sex + sch_type*ethnic + raven + reading + maths + lem_total,
  reading = reading ~ sex + sch_type*ethnic + raven + b_total + maths + lem_total,
  maths = maths ~ sex + sch_type*ethnic + raven + b_total + reading + lem_total,
  lem_total = lem_total ~ sex + sch_type*ethnic + raven + b_total + reading + maths
)

# Imputace pomocí predictive mean matching
df_imp_pmm <- mice(df_nabular, 
                   method = "pmm", # predictive mean matching imputace
                   formulas = form_list, # Specifikace imputačních modelů.
                   auxilary = FALSE, # nepřidávat jako prediktory žádné další proměnné
                   donors = 5, # počet potenciálních donors
                   m = 10, # Počet multiple imputations 
                   maxit = 10, # Počet iterací pro "ustálení" imputovaných hodnot
                   seed = 123) # Iniciální hodnota generátoru náhodných čísel

# Takto sloučíme všechny imputované datasety do jednoho v longer formátu. 
df_imp_pmm_long <- df_imp_pmm %>% 
  complete(action = "long") %>% 
  as_tibble()
df_imp_pmm_long

# Můžeme si dále vytvoři sloupce indikující, zda v lem_total nebo reading 
# byla původně chybějící hodnota.
df_imp_pmm_long <- df_imp_pmm_long %>% 
  mutate(lem_or_reading_na = lem_total_NA == "NA" | reading_NA == "NA")

df_imp_pmm_long %>% 
  ggplot(aes(lem_total, reading, color = lem_or_reading_na)) +
  geom_point() +
  facet_wrap(~.imp, ncol = 2) +
  theme(legend.position = "top")

# Otázka: Poznali byste nyní (kdyby nebylo barevného odlišení), 
# které hodnoty byly imputovány a které nikoli?

# Také bychom měli ověřit, že se imputované hodnoty “ustálily”. 
# Můžeme to ověřit nepřímo např. tak, že se podíváme, 
# zda průměry a směrodatné odchylky imputovaných proměnných jsou stabilní 
plot(df_imp_pmm)  

# Nakonec si můžeme ukázat, jak imputované datasety použít k odhadu 
# lineárně regresního modelu. 
fit <- with(
  df_imp_pmm,
  lm(reading ~ sex + sch_type*ethnic + lem_total)
)

# term = prediktor;
# m = počet imputací;
# estimate = odhad koeficientu modelu;
# ubar = průměrný vnitroimputační rozptyl (rozptylem se zde myslí error 
#        variance čili čtverec standardní chyby);
# b = meziimputační rozptyl;
# t = celkový rozptyl souhrného odhadu (ubar + b) čili 
#     celkový čtverec standardní chyby;
# dfcom = reziduální stupně volnosti, pokud bychom měli kompletní data;
# df = reziduální stupně volnosti (pro testy regresních koeficientů);
# riv = odhad relativního nárůstu chybového rozptylu v důsledku chybějících dat;
# lambda = Odhad podílu chybového rozptylu způsobeného chybějícími daty;
# fmi = Rovněž odhad podílu chybového rozptylu způsobeného chybějícími daty, 
#       ale s korekcí podhodnocení (které je dáno tím, že namáme nekonečně 
#       mnoho imputovaných datasetů)


# Odhady standardních chyb, testové statistiky a p-hodnot pak dostaneme, 
# pokud output funkce pool() pošleme do funkce summary().
summary(pool(fit), conf.int = TRUE, conf.level = .95) %>% 
  mutate(across(where(is.double), ~round(.x, digits = 3)))

pool.r.squared(fit)
pool.r.squared(fit_pmm2, adjusted = TRUE)

# Kdybychom měli více nested modelů, můžeme je srovnat pomocí funkce
# mice::D3(), 
# která provádí likelihood-ratio test.