##%######################################################%##
#                                                          #
####                    Missing data                    ####
#                                                          #
##%######################################################%##


# Preparation -------------------------------------------------------------
library(mice)
library(tidyverse)
library(naniar)
library(simputation)
library(VIM)
source("R101_funs.R")

# Create x, a vector, with values NA, NaN, Inf, ".", and "missing"
x <- c(NA, NaN, Inf, ".", "missing")

# Jednoduché funkce pro zjištění výskytu chybějících hodnot
any_na(x) # Je jakákoliv hodnota NA?
are_na(x) # Které prvky jsou NA?
is.na(x)

# Použitý dataset
df <- foreign::read.spss("data/lem.sav", to.data.frame = TRUE,
                         trim.factor.names = TRUE, reencode = TRUE) %>% 
  as_tibble() %>% 
  mutate(
    across(c(school, class, edu_a), as.factor),
    b_total = b1 + b2
  )


# Přehled použitého datasetu
glimpse(df)
summary(df)
attributes(df)$variable.labels

# n_miss() zjistí celkový počet chybějících hodnot napříč proměnnými
n_miss(df)

# funkce prop_miss() a prop_complete() zjistí podíl chybějících a nechybějících hodnot
prop_miss(df)
prop_complete(df)

# Tabulky -----------------------------------------------------------------
# Funkce miss_var_summary() zjistí počet a podíl chybějících hodnot 
# v jednotlivých proměnných
# A funkce miss_var_table() ukáže kolika proměnným chybí kolik hodnot
miss_var_summary(df)
miss_var_table(df)

# Pomocí grou by můžeme srovnat různé skupiny
df %>% 
  group_by(ethnic) %>% 
  miss_var_summary() %>% 
  pivot_wider(names_from = ethnic, values_from = c(n_miss, pct_miss))

# Funkce miss_case_summary() ukáže počet NAs pro každého respondenta
# a funkce miss_case_table() kolika respondentům chybí kolik hodnot
miss_case_summary(df)
miss_case_table(df)

df %>% 
  slice(12:17)

# Funkcí add_n_miss() můžete přidat do datasetu sloupec s počtem chybějících 
# hodnot v daném řádku
df %>% 
  add_n_miss() %>% 
  select(n_miss_all, everything())


miss_by_group(df, ethnic)
miss_by_group(df, sex, ethnic)

# Grafy -------------------------------------------------------------------
# Která proměnná má nejvíce chybějících hodnot?
gg_miss_var(df)
gg_miss_var(df, facet = ethnic, show_pct = TRUE)
gg_miss_var(df, facet = sex, show_pct = TRUE)

# Shluky proměnných, jejichž hodnoty často chybí společně
vis_miss(df, cluster = TRUE) + 
  theme(axis.text.x = element_text(angle = 90))


# Nejčetnější vzorce chybějících hodnot
# Let’s find out.
gg_miss_upset(df)
gg_miss_upset(df, 
              nsets = 10, # max počet proměnných
              nintersects = 10) # počet vzorců


# Heatmapa chybějících hodnot pro různé skupiny
gg_miss_fct(df, fct = ethnic)
gg_miss_fct(df, fct = school)
gg_miss_fct(df, fct = sch_type)

gg_miss_fct(df, fct = .data$sch_type)

# Scatterploty s chybějícími hodnotami
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))


# Nabular data ------------------------------------------------------------
# Shadow matrix (stínová matice) `as_shadow()`
# Obsahuje jednoduše údaje o tom, zda je hodnota chybějící (NA) nebo nechybějící (!NA)
as_shadow(df) 

#  bind_shadow() vytvoří shadow matrix a připojí ji jako nové sloupce k původním datům
bind_shadow(df)

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

# Shadow matrix pak můžeme využít porovnání skupin vymezených podle 
# chybějících dat v nějaké proměnné
df_shadow %>%
  group_by(raven_NA) %>% 
  summarize(
    n = n(),
    M_lem = mean(lem_total, na.rm = TRUE),
    SD_lem = sd(lem_total, na.rm = TRUE))

df_shadow %>%
  group_by(maths_NA) %>% 
  summarize(
    n = n(),
    M_lem = mean(lem_total, na.rm = TRUE),
    SD_lem = sd(lem_total, na.rm = TRUE))

df_shadow %>%
  group_by(edu_f_NA) %>% 
  summarize(
    n = n(),
    M_lem = mean(lem_total, na.rm = TRUE),
    SD_lem = sd(lem_total, na.rm = TRUE))

# Graf hustoty pravděpodobnosti podle chybějí v reading a podle etnika
df_shadow %>% 
  ggplot(aes(lem_total, color = reading_NA)) +
  geom_density() +
  facet_wrap(~ethnic)

# 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 tedy, že v základní škole 
# se dítě naučí více než ve zvláštní. V duchu projektu budeme zvažovat, 
# že efekt š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ě velmi ovlivněna.


df_shadow <- df %>% 
  bind_shadow() %>% # připojení shadow matrix
  impute_below_all() %>% # imputace chybějících hodnot mimo rozmezí pozorovaných hodnot
  add_label_shadow() # přidat proměnnou "any_missing"

ggplot(df_shadow, aes(lem_total, reading, color = any_missing)) +
  geom_point()

# Doporučení pro tvorbu imputačního modelu
# 1. Zahrňme všechny proměnné které budou součastní finálního modelu, 
# jejž chceme fitovat a to včetně závislé proměnné (ta bude také použita 
# jako predikor chybějících dat v ostatních proměnných) a odůvodněných
# interakcí. Jinak mohou být odhady finálního modelu zkreslené 
# (směrem k nule, tedy podhodnocené)

# 2. Zahrňe všechny proměnné, které souvisejí s absencí odpovědí v cílových
# proměnných. Můžeme tak učinit na základě teoretické úvahy a na empirickém 
# základě. To znamená, že bychom měli zahrnout proměnné, jejichž distribuce 
# se výrazně liší mezi skupinami s/bez chybějících dat. Obvykle to můžeme 
# zjistit, když blíže prozkoumáme korelace různých proměnných s indikátorem
# chybějících hodnot (!NA vs NA) v cílovových (imputovaných) proměnné.

# 3. Zahrňme proměnné vysvětlující významné podíl rozptylu v proměnných,
# které plánujem imputovat. Tím snižujeme variabilitu imputovaných hodnot
# I v tomto případě můžeme prozkoumat, jak korelují ostatní proměnné
# s cílovými (imputovanými) proměnnými.

# 4. Vyřaďte ty proměnné z kroku 2 a 3, které mají příliš chybějících hodnot 
# v podskupině nekompletních případů (tj. případů, které mají nějakou
# chybějící hodnotu v proměnné cílového modelu)
# Jednoduchým ukazatelem je podíl chybějících hodnot v této skupině.

# * Mean imputation -------------------------------------------------------

# Jednoduchá imputace průměru
df_imp_mean <- df %>% 
  bind_shadow() %>% 
  add_label_shadow() %>% 
  mutate(
    across(where(is.numeric), impute_mean) # imputovat průměr do všech numerických proměnných
  )

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
df %>% 
  select(where(is.numeric)) %>% 
  cor(use = "pairwise.complete.obs") %>% 
  round(digits = 2) %>% 
  .[, c("reading", "lem_total")]

df_imp_mean %>% 
  select(where(is.numeric)) %>% 
  cor(use = "pairwise.complete.obs") %>% 
  round(digits = 2)%>% 
  .[, c("reading", "lem_total")]


# * Hot deck imputation -----------------------------------------------------
miss_var_summary(df) %>% view()

df_imp_hd <- df %>% 
  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 = "maths") # seřadit podle

# Draw a margin plot
df_imp_hd %>% 
  select(lem_total, reading, ends_with("_imp")) %>% 
  marginplot(delimiter = "_imp")



# * kNN imputation ----------------------------------------------------------
df_imp_knn <- df %>% 
  kNN(variable = c("lem_total", "reading"),
      dist_var = c("sch_type", "maths", "raven",
                   "b1", "b2"),
      k = 5, # počet neighbors
      numFun = weighted.mean,
      weightDist = TRUE)

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


# * Regression imputation -------------------------------------------------
# Funkce pro výpočet průměrné absolutní procentuální změny (mezi iteracemi)
# MAPC = mean absolute percent change
mapc <- function(a, b) {
  mean(abs(b - a) / a, na.rm = TRUE)
}

diff_lem_total <- double(10)
diff_reading <- double(10)

na_lem_total <- is.na(df$lem_total)
na_reading <- is.na(df$reading)

# Inicializace chybějících hodnot pomocí prosté hotdeck imputece
df_imp_lm <- df %>% 
  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 = "maths") # seřadit podle

for (i in 1:10) {
  # Uložení předchozí iterace
  prev_iter <- df_imp_lm
  # Nastavit na NA hodnoty, které byly původně NA
  df_imp_lm$reading[na_reading] <- NA
  # Imputovat reading pomocí regrese (predikované hodnoty)
  df_imp_lm <- df_imp_lm %>% 
    impute_lm(reading ~ sch_type + ethnic + lem_total + sex)
  
  # Nastavit na NA hodnoty, které byly původně NA
  df_imp_lm$lem_total[na_lem_total] <- NA
  # Imputovat lem_total pomocí regrese (predikované hodnoty)
  df_imp_lm <- df_imp_lm %>% 
    impute_lm(lem_total ~ sch_type + ethnic + reading + sex)
  # Výpočet MAPC
  diff_reading[i] <- mapc(prev_iter$reading, df_imp_lm$reading)
  diff_lem_total[i] <- mapc(prev_iter$lem_total, df_imp_lm$lem_total)
}

plot(diff_reading, type = "b")
plot(diff_lem_total, type = "b")

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

# Multiple imputation with mice -------------------------------------------
df_shadow <- df %>% 
  bind_shadow() %>% 
  add_any_miss(lem_total, reading, sex, sch_type, ethnic)

df_shadow

fit_listwise <- lm(reading ~ sex + sch_type*ethnic + lem_total, data = df)
summary(fit_listwise)


# Opět pomocí lineární regrese.
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
)

# Impute the missing data with regression imputation
df_imp_lm <- mice(df_shadow, 
                  method = "norm.predict", # linear regression, predicted values
                  formulas = form_list, 
                  auxiliary = FALSE, # Nepřidávat jako prediktory žádné další proměnné
                  m = 1, # Počet multiple imputations 
                  maxit = 5) # Počet iterací
df_imp_lm

complete(df_imp_lm) %>% 
  ggplot(aes(lem_total, reading, color = any_miss_vars)) +
  geom_point()

fit_lm_pv <- lm(reading ~ sex + sch_type*ethnic + lem_total, data = complete(df_imp_lm))

summary(fit_listwise)
summary(fit_lm_pv)


# * Predictive mean matching ------------------------------------------------
df_imp_pmm <- mice(df_shadow, 
                   method = "pmm", # linear regression, predicted values
                   formulas = form_list, 
                   auxiliary = 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 = 5, # Počet iterací
                   seed = 123)

# Scatterploty pro všechny imputace
c(1:10) %>% 
  map(~complete(df_imp_pmm, action = .x)) %>% 
  map(~ggplot(.x, aes(lem_total, reading, color = any_miss_vars)) +
       geom_point())

# Grafy hustoty pravděpodobnosti pro všechny imputace
c(1:10) %>% 
  map(~complete(df_imp_pmm, action = .x)) %>% 
  map(~ggplot(.x, aes(lem_total, color = any_miss_vars)) +
        geom_density())

c(1:10) %>% 
  map(~complete(df_imp_pmm, action = .x)) %>% 
  map(~ggplot(.x, aes(reading, color = any_miss_vars)) +
        geom_density())

# Fit linear regression to each imputed data set
# Model 1
fit_pmm1 <- with(
  df_imp_pmm,
  lm(reading ~ sex + sch_type*ethnic)
)
# Model 2
fit_pmm2 <- with(
  df_imp_pmm,
  lm(reading ~ sex + sch_type*ethnic + lem_total)
)

# Pool regression results
pool(fit_pmm1)
pool(fit_pmm2)
# ubar = The mean of the variances (i.e. the pooled within-imputation variance)
# b = The between-imputation variance
# t = The total variance of the pooled estimated
# The relative increase in variance due to nonresponse
# df = Residual degrees of freedom for hypothesis testing.
# riv = Relative increase in variance due to nonresponse
# lambda = Proportion of total variance due to missingness
# fmi = Fraction of missing information
pool(fit_pmm1) %>% summary() 
pool(fit_pmm2) %>% summary() 

pool.r.squared(fit_pmm1)
pool.r.squared(fit_pmm2)

summary(pool(fit_pmm1), conf.int = TRUE, conf.level = .995)
summary(pool(fit_pmm2), conf.int = TRUE, conf.level = .995)

# LR test pro srovnání modelů
D3(fit_pmm2, fit_pmm1)