# Generisemo podatke iz visedimenzione normalne raspodele,
# takve da obe kolone imaju srednju vrednost 3 i standardnu
# devijaciju 1. Takodje, postavljamo kovarijaciju na 0.6, kako bismo imali
# zavisnost izmedju njih,
set.seed(216)
n <- 1000
cov_mat <- matrix(c(1, 0.6, 0.6, 1), nrow = 2)
library(MASS)
X <- mvrnorm(n, mu = c(3, 3), Sigma = cov_mat)
plot(X) #postoji zavisnost izmedju X1 i X2

x1 <- X[,1]
x2 <- X[,2]
mean(x2)
sd(x2)

# Videcemo kako pojava nedostajucih podataka moze na razlicite nacina da utice
# na prirodu podataka. Neodredjena polja (NA) ubacujemo u x2, i posebno cemo
# obratiti paznju na koji nacin to utice na vrednosti uzoracke sredine i
# standardne devijacije.

# MCAR - potpuno nasumicno nedostajanje
# Ubacujemo nedostajuce vrednosti u podatke potpuno nasumicno, u svakom redu
# postoji 50% sansa da podataka nestane
missing <- sample(c(TRUE, FALSE), n, replace = TRUE)
x2_mcar <- x2
x2_mcar[missing] = NA
mean(x2_mcar) 
sd(x2_mcar)
# Ne mogu da se izracunaju, izbacicemo nedostajuce vrednosti
x2_mcar_obs = na.omit(x2_mcar)
mean(x2_mcar_obs)
sd(x2_mcar_obs)
par(mfrow = c(2,1))
boxplot(x2_mcar_obs, ylim = c(0,7), horizontal = TRUE, main = 'Podaci sa nedostajucim vrednostima')
boxplot(x2, ylim = c(0,7), horizontal = TRUE, main = 'Originalni podaci')
# Ocene ostaju nepristrasne, samo su manje precizne jer imamo samo 60%
# originalnih podataka

# MAR - nasumicno nedostajanje
# U ovom slucaju dozvoljavamo da vrednost x1 utice na to da li ce podatak za x2
# nedostajati

# Probacemo sa raznim varijantama

# Desno nedostajanje
# p(X2 = NA|X1) = sigma(x1 - 4), gde je sigma logisticka funkcija
# Vece vrednosti x1 daju vecu verovatnocu nedostajanja za x2
# Posto su x1 i x2 pozitivno korelisani, ovo znaci da ce cesce da nedostaju
# podaci za x2 kad je x2 veci, dajuci prirodno pristrasnost nasim ocenama
prob_of_missing <- plogis(x1-3)
missing <- as.logical(rbinom(n, 1, prob = prob_of_missing))
x2_mar_right = x2
x2_mar_right[missing] = NA
x2_mar_right_obs = na.omit(x2_mar_right)
mean(x2_mar_right_obs)
sd(x2_mar_right_obs)
boxplot(x2_mar_right_obs, ylim = c(0,7), horizontal = TRUE, main = 'Podaci sa nedostajucim vrednostima')
boxplot(x2, ylim = c(0,7), horizontal = TRUE, main = 'Originalni podaci')
# Sada su obzervabilni podaci pristrasno ulevo

# Centralno nedostajanje:
prob_of_missing <- plogis(0.75 - abs(x1 - 3))
missing <- as.logical(rbinom(n, 1, prob = prob_of_missing))
x2_mar_centre = x2
x2_mar_centre[missing] = NA
x2_mar_centre_obs = na.omit(x2_mar_centre)
mean(x2_mar_centre_obs)
sd(x2_mar_centre_obs)
boxplot(x2_mar_centre_obs, ylim = c(0,7), horizontal = TRUE, main = 'Podaci sa nedostajucim vrednostima')
boxplot(x2, ylim = c(0,7), horizontal = TRUE, main = 'Originalni podaci')
# Posto je verovatnoca izbacivanja bila veca kad je x1 blize centru, sada su
# nam podaci vestacki rasprseniji, standardna devijacija je pozitivno pristrasna

# Nedostajanje na repovima:
prob_of_missing <- plogis(abs(x1 - 3))
missing <- as.logical(rbinom(n, 1, prob = prob_of_missing))
x2_mar_tail = x2
x2_mar_tail[missing] = NA
x2_mar_tail_obs = na.omit(x2_mar_tail)
mean(x2_mar_tail_obs)
sd(x2_mar_tail_obs)
boxplot(x2_mar_tail_obs, ylim = c(0,7), horizontal = TRUE, main = 'Podaci sa nedostajucim vrednostima')
boxplot(x2, ylim = c(0,7), horizontal = TRUE, main = 'Originalni podaci')
# Posto je ovde verovatnoca izbacivanja bila veca za tacke na repu, sada su
# podaci vestacki suzeni, standardna devijacija je negativno pristrasna

# MNAR - Nenasumicno nedostajanje
# U ovom slucaju vrednost x2 takodje utice na svoje nedostajanje - ovo je
# sustinski neotklonjiv problem

# Desno nenasumicno nedostajanje:
prob_of_missing <- plogis(x2-3)
missing <- as.logical(rbinom(n, 1, prob = prob_of_missing))
x2_mnar_right = x2
x2_mnar_right[missing] = NA
x2_mnar_right_obs = na.omit(x2_mnar_right)
mean(x2_mnar_right_obs)
sd(x2_mnar_right_obs)
boxplot(x2_mnar_right_obs, ylim = c(0,7), horizontal = TRUE, main = 'Podaci sa nedostajucim vrednostima')
boxplot(x2, ylim = c(0,7), horizontal = TRUE, main = 'Originalni podaci')
# Podaci su pomereni ulevo, ali jos je naglasenije

# U svim ovim slucajevima, osim sto smo ubacili pristrasnost u podatke, bilo
# u standardnoj devijaciji ili srednjoj vrednosti, te ocene su takodje bile
# manje precizne, odnosno, imale su vecu disperziju, zbog cinjenice da smo
# izbacivali znacajan deo podataka.



## Imputacija nedostajucih vrednosti
# Umesto izbacivanja podataka, probacemo sa metodama imputacije podataka
# Probacemo nedostajuce vrednosti da zamenimo necim smislenim

# Imputacija srednjom vrednoscu/medijanom
# Radicemo na primeru desnog nedostajanja
mean(x2_mar_right_obs)
missing = is.na(x2_mar_right)

x2_mar_right_imputed_mean = x2_mar_right
x2_mar_right_imputed_mean[missing] = mean(x2_mar_right_obs)
mean(x2_mar_right_imputed_mean)
sd(x2_mar_right_imputed_mean)
# Srednja vrednost podataka naravno ostaje ista (i dalje pristrasna), ali
# smo takodje smanjili varijabilnost u podacima time sto smo sve nedostajuce
# podatke zamenili srednjom vrednoscu.
# Jedina prednost ovoga u odnosu na izbacivanje podataka je da ako, npr.
# x1 i x2 predstavljaju prediktore nekog modela, ne izbacujemo celu opservaciju
# samo zato sto ne znamo koja je vrednost x2.

# Imputacija regresijom
# Pravimo model x2 ~ x1, i time mozemo u MAR slucaju bolje da ocenimo x2
# Model pravimo nad podacima koji su nam dostupni
X = data.frame(x1 = x1,x2 = x2_mar_right)
X_obs = X[!missing,]
X_miss = X[missing,]
model <- lm(x2 ~ x1, data = X_obs)
summary(model) 
# R^2 = 0.27, nije mnogo jaka veza, ali nam x1 svakako daje dodatnu informaciju
# o vrednosti x2
x2_mar_right_imputed_lm = x2_mar_right
x2_mar_right_imputed_lm[missing] = predict(model, newdata = X_miss)
mean(x2_mar_right_imputed_lm)
sd(x2_mar_right_imputed_lm)
# Sada x2 vise nije pomeren ulevo, resili smo problem pristrasnosti srednje
# vrednosti. Ali opet imamo problem da smo vestacki smanjili varijabilnost
# podataka.


# Predictive Mean Matching (PMM) i visestruka imputacija
# Predictive Mean Matching isto pravi linearni model na osnovu ostalih ostalih
# prediktora, ali ima odredjenu nasumicnost pri odabiru vrednosti. Osim toga,
# robusniji je metod od linearne regresije, bolje radi kada veza nije bas
# linearna.
# Zbog nasumicnosti ove metode, praksa je da se radi visestruka imputacija:
# Umesto da NA vrednosti menjamo jedinstvenom vrednoscu, pravimo vise baza
# od nase pocetne, gde svaka ima razlicite imputirane vrednosti
x2_mar_right[missing] = NA

library(mice)
# Funkcija mice sluzi za visestruku imputaciju, a pretpostavljeni metod je
# PMM, broj imputacije je pretpostavljeno 5
imputed_data = mice(X, m = 5, method = 'pmm')
imputed_datasets = complete(imputed_data, 'all') #vraca listu baza
# Mozemo da vidimo sta smo dobili, originalna baza ima NA vrednosti u koloni x2,
# a mi sad imamo 5 baza koje su popunili te vrednosti na razlicite nacine
View(X)
View(imputed_datasets[[1]])
View(imputed_datasets[[2]])

# Prikazujemo srednje vrednosti i standardne devijacije za svih 5 baza i
# vidimo da su priblizno jednake originalnim podacima bez nedostajanja.
means <- sapply(imputed_datasets, function(df) mean(df[,2]))
means
sds <- sapply(imputed_datasets, function(df) sd(df[,2]))
sds

# Mozemo i sa grafika da vidimo kako PMM zadrzava strukturu podataka:
# Nacrtacemo kako izgledaju podaci sa razlicitim nacinima imputacije
# Crnom su obojeni opservabilni podaci, crvenom imputirani
par(mfrow = c(1,1))
plot(x1, x2_mar_right_imputed_mean, col = missing + 1, pch = 16, main = "Imputacija srednjom vrednoscu")
plot(x1, x2_mar_right_imputed_lm, col = missing + 1, pch = 16, main = "Imputacija linearnom regresijom")
plot(imputed_datasets[[1]], col = missing + 1, pch = 16, main = "Imputacija sa Predictive Mean Matching")


## Visestruka imputacija pri pravljenju modela
library(VIM)
data(sleep)
?sleep
View(sleep)

nrow(sleep)
nrow(na.omit(sleep))
# Izbacivanjem podataka bismo izgubili 20 redova, sto je skoro trecina

# Obrazac nedostajucih podataka mozemo da prikazemo na sledeci nacin:
md.pattern(sleep)
# U 42 reda imamo sve vrednosti
# U 9 imamo nedostajanja u Dream i NonDream
# U 3 imamo samo nedostajanje u Gest,itd.

# Posto cemo da predvidjamo kolicinu sna, nema smisla da imamo podatke o sanjanju
data = sleep[,-c(3,4)]

# Posto predvidjamo Sleep, svakako zelimo da izbacimo redove gde nemamo taj
# podatak
data = data[!is.na(sleep$Sleep),]

# Posto model treba da predvidja Sleep, ne sme te podatke da koristi u procesu
# imputacije
X = data[,-3]
y = data[,3]
imp = mice(X, m = 5) #5 imputacija
X_imputed = complete(imp, "all")
# Sada u X_imputed imamo 5 baza, nad svakom pravimo model
models <- lapply(X_imputed, function(df) lm(y ~., data = df))

# U mice paketu imamo i funkciju pool, koja ove modele kombinuje
final_model <- pool(models)
summary(final_model)
# Vidimo koeficijente modela, kao i p-vrednosti. Nije bas lepo kao summary
# za lm, nemamo zvezdice, nemamo R^2

# Pravilo za dobijanje koeficijenata finalnog modela je jednostavno, samo se
# uzme prosek koeficijenata svakog modela
coefs = sapply(models, function(model) model$coefficients)
#coefs
final_coefs = rowMeans(coefs) 
final_coefs #isto kao u summary(final_model)

# Sada cemo sami da izracunamo R^2 ovog modela
# Predikciju takodje radimo nad imputiranim podacima, i uzimamo prosek
# Imamo i problem da predict funkcija ne radi sa finalnim modelom, pa cemo
# rucno racunati predikcije, koristeci koeficijente

predict_with_coefs <- function(coefs, newdata)
  apply(newdata, 1, function(x) {
    coefs[1] + sum(coefs[-1]*x) #slobodan clan + skalarni proizvod ostalih koeficijenata i vrste u podacima
  })

y_pred = sapply(X_imputed, function(df) predict_with_coefs(final_coefs, df))
head(y_pred) #Predikcija nad svakom bazom
y_pred_final = rowMeans(y_pred)

# Srednjekvadratna greska
MSE = mean((y - y_pred_final)^2)
n <- length(y)
R2 = 1 - (MSE*n)/(var(y)*(n-1))
MSE
R2


# Navedeni proces visestruke imputacije je jako bitan kada nam je u interesu
# da donosimo statisticke zakljucke nad bazama podataka - npr. ako zelimo da
# vrsimo selekciju prediktora na osnovu p-vrednosti testova znacajnosti.
# Ipak, ovaj proces je pipav (funkcija pool ne radi za bilo koji model), i
# cesto naporan (ponavljamo transformacije, analize, itd. nad svakom imputiranom
# bazom).

# Moj predlog za one koje ne zele da prihvate izazov visestruke imputacije je
# sledeci:
# Uzmemo medijane svih imputiranih vrednosti i spojimo ih u jednu bazu. Moramo
# da budemo svesni samo da gubimo ocuvanost strukture podataka ovim. Konkretno
# vestacki smanjujemo varijaciju u obelezjima sa nedostajanjima.
# Kada ovo radimo, veca vrednost m (vise imputiranih vrednosti) ce dati
# precizniju imputaciju, ali ce losije odrzavati originalnu strukturu podataka.
# Ova metologija nece dati jednako dobre modele na kraju, ali se u praksi
# cesce koristi jer znacajno jednostavnija.
pmm_median = function(X, m = 5)
{
  imputation = mice(X, m = m, method = 'pmm')
  X_imputed = complete(imputation, 'all')
  X_imputed_pooled = X_imputed[[1]] # za pocetak uzimamo prvu bazu
  for(i in 1:nrow(X_imputed_pooled))
    for(j in 1:ncol(X_imputed_pooled))
    {
      # Za svaku bazu izvlacimo vrednost i-te vrste i j-te kolone:
      sve_vrednosti = c()
      for(l in 1:m)
        sve_vrednosti[l] = X_imputed[[l]][i,j]
      # Uzimamo medijanu:
      X_imputed_pooled[i,j] = median(sve_vrednosti)
    }
  return(X_imputed_pooled)
}

# Isprobajmo ovo na generisanim podacima od ranije:
X_sim = data.frame(x1 = x1,x2 = x2_mar_right)
X_sim_imp = pmm_median(X_sim)
mean(X_sim_imp[,2])
sd(X_sim_imp[,2]) #manje od originalnih ~1.06
plot(X_sim_imp, col = missing + 1, pch = 16, main = "Imputacija sa medijanom Predictive Mean Matching")
# Crvene imputirane vrednosti su malo manje rasprsene od originalnih crnih.

# Vratimo se na model spavanja sisara:
X_imp_pooled = pmm_median(X)
# Dalje radimo kao sa obicnom bazom, s tim sto moramo da pazimo da p-vrednosti
# testova nisu bas verodostojne.
model = lm(y ~., data = X_imp_pooled)
summary(model)