# Generisemo podatke iz visedimenzione normalne raspodele,
# takve da obe kolone imaju srednju vrednost 3 i standardnu
# devijaciju 1. Takod
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
# 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]])

means <- sapply(imputed_datasets, function(df) mean(df[,2]))
means
sds <- sapply(imputed_datasets, function(df) sd(df[,2]))
sds
# Vidimo, na osnovu srednjih vrednosti, da podaci vise nisu pomereni ulevo, a
# pritom smo zadrzali varijabilnost u podacima.
# Posto smo sada imamo 5 baza, finalne ocene za srednju vrednost i standardnu
# devijaciju mozemo da dobijemo uprosecavanjem ocena iz baza.
mean_final = mean(means)
sd_final = mean(sds)
mean_final
sd_final

# Mozemo i sa grafika da vidimo kako PMM zadrzava strukturu podataka:
# Nacrtacemo kako izgledaju podaci sa razlicitim nacinima imputacije
# Crnom su obojeni obzervabilni 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")


## 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

# Posto cemo da predvidjamo kolicinu sna, nema smisla da imamo podatke o snovima
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)
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