## Logisticka regresija
# Sada cemo se baviti logistickom regresijom sa vise od jednog prediktora

# Koristicemo podatke o primanjima amerikanaca, pri cemu cemo pokusati da
# modelujemo da li odredjena osoba ima primanja preko $50K godisnje na osnovu
# drugih faktora.
# Podaci: http://www.matf.bg.ac.rs/p/files/94-income.csv
income <- read.csv("income.csv")
# plot(income)
# Napravicemo model koji modelira ABOVE50K preko ostalih promenljivih
logistic <- glm(ABOVE50K ~ ., income, family = binomial)
summary(logistic)
# Statisticka znacajnost pojedinacnih koeficijenata se testira Valdovim testom,
# koji se oslanja na normalnost raspodele ocene maksicmalne verodostojnosti. Iz
# pregleda modela, izgleda da su svi koeficijenti znacajni, osim RACE. To cemo
# analizirati dodatno.

# Pogledajmo prvo neke vrednosti koje se dobijaju. 
# Vise nemamo R^2 kao ocenu kvaliteta, vec nam pise vrednost AIC(Akaike
# Information Criterion), koji mozemo koristiti za poredjenje razlicitih modela
# - sto je manji AIC to je bolji model. Ne postoji "referentna" AIC vrednost
# koja se moze koristiti za merenje kvaliteta jednog modela, vec se koristi kao
# mera za poredjenje dva modela.
# AIC se definise kao -2*logLikelihood(model) + 2*(p+1)
logL <- as.numeric(logLik(logistic))
-2*logL + 2*length(logistic$coefficients)
# Za nas model AIC je oko 5067

# Takodje imamo vrednosti za Null deviance i Residual deviance (pogledati
# teoriju za formule). Ove vrednosti su znacajne za testiranje da li postoji
# ikakav uticaj prediktora na zavisnu promenljivu, poput F testa u linearnoj
# regresiji. Ako su Null i Residual deviance slicne, to ukazuje na to da ne
# postoji znacajan uticaj prediktora na zavisnu promenljivu.

# Deviance se racuna kao razlika logaritma verodostojnosti zasicenog modela i
# modela koji posmatramo. Zasicen model je onaj model koji ima isti broj
# prediktora kao i podataka, tj. u optpunosti odgovara podacima.

# Null deviance je devijance za model samo sa slobodnim clanom i racuna se kao:
null_model <- glm(ABOVE50K ~ 1, income, family = binomial)
-2*logLik(null_model)
# Residual deviance je deviance za nas model
-2*logLik(logistic)

# Testiranje znacajnosti razlika izmedju deviance, tj. znacajnosti razlika
# izmedju dva (ugnjezdjena) modela se vrsi pomocu anova funkcije, kao kod
# linearnih modela, s tim sto se prosledjuje argument test="Chisq". Ovo
# testiranje predstavlja test kolicnika verodostojnosti. (videti teoriju)

# Testiranje da li postoji zncajan uticaj prediktora na zavisnu promenljivu,
# poput F testa od pre, se vrsi sa
anova(null_model, logistic, test = "Chisq")
# Znacajna je vrednost test statistike, pa smatramo da nisu svi koeficijenti 0.

# Uvideli smo ranije da znacaj prediktora vezanih za rasu nije veliki. Zato cemo
# testirati koristeci test kolicnika verodostojnosti znacajnost prediktora RACE.
no_race <- update(logistic, . ~ - RACE) # menjamo logistic da bude bez RACE
anova(no_race, logistic, test = "Chisq")
# Ipak je p vrednost vrlo mala! Postoji bitan uticaj rase na primanja, iako su
# svi koeficijenti bili beznacajni.

# Primetimo da je ovde referentni nivo rase bio Amer-Indian-Eskimo,
income$RACE
# a bolje bi bilo da bude White, jer je to vecinska populacija. Zato cemo to
# promeniti i videti sta se desava.
income$RACE <- relevel(income$RACE, " White")

# Novi model postaje
logistic <- glm(ABOVE50K ~ ., income, family = binomial)
summary(logistic)
# Sada je koeficijent uz RACE Black vrlo znacajan i pritom je negativan, sto
# znaci da crnci imaju znacajno manje sanse da imaju primanja veca od $50K
# godisnje. Pitanje je da li je ovo posledica diskriminacije ili nekog drugog
# problema u drustvu, ali time se nebavimo. Ovo je bio razlog zasto nam test
# kolicnika verodostojnosti govori da postoji znacjan uticaj rase na model.

no_race <- update(logistic, . ~ - RACE) # menjamo logistic da bude bez RACE
anova(no_race, logistic, test = "Chisq")
# Vrednost ovog testa se nije promenila od ranije, vec samo testovi za
# pojedinacne koeficijente.

# Interesantna velicina su deviance reziduali, jer je suma njihovih kvadrata
# jednaka upravo reisudal deviance i metodom maksimalne verodostojnosti mi
# minimizujemo upravo ovu sumu kvadrata deviance reziduala. U kontekstu
# masinskog ucenje, ovi reziduali se koriste kao fukncija gubitka za logisticku
# regresiju. Dobijaju se ovako...
dev <- residuals(logistic, type="deviance")
# Anjihova suma kvadrata je bas jednaka residual deviance
sum(dev^2)

# Za zadatak klasifikacije, mozemo odrediti povoljan prag za verovatnocu pomocu
# skupa za testiranje datog u income_test.csv. Trudicemo se da nadjemo prag koji
# maksimizuje F1 skor.
# Podaci: http://www.matf.bg.ac.rs/p/files/94-income_test.csv
library(Metrics)
income_test <- read.csv("income_test.csv")
thresholds <- seq(0.1, 0.9, 0.1)
scores <- sapply(thresholds, function(thr){
  preds <- predict(logistic, income_test, type="response") > thr
  fbeta_score(income_test$ABOVE50K, preds)
})
scores
thresholds[which.max(scores)]
# Najbolje ponasanje imamo za prag 0.4

library(pROC)
plot(roc(income_test$ABOVE50K,
     predict(logistic, income_test, type="response")),
     print.auc = TRUE)
# AUC nam je 0.82 sto je vrlo dobro