# Logisticka regresija

# Teorijski uvod je u dodatku na sajtu.

# Koristicemo iris skup podataka.
# nacrtajmo sta sve ima u njemu
plot(iris)
# Izvucicemo samo 2 vrste cvetova
iris <- iris[iris$Species != "setosa", ]
iris$Species <- factor(iris$Species)
plot(iris)
# Nas zadatak ce biti da klasifikujemo cvetove po vrsti, koristeci ostale
# promenljive u podacima.

# Koristimo logisticku regresiju da modelujemo verovatnocu da je cvet vrste
# "virginica" (iz contrasts(iris$Species vidimo da je sa 1 oznacena ova vrsta)

# Uopsteni linearni modeli se prave funkcijom glm, gde je najbitniji parametar
# "family" koji oznacava koju familiju raspodela ima zavisna promenljiva. Kod
# logisticke regresije se koristi family=binomial jer pretpostavljamo binomnu
# raspodelu za Y (bernulijeva je binomna). Koristicemo kanonsku link funkcju
# (default vrednost u R) "logit"

# Uzecemo u obzir samo prediktor Petal.Width
model <- glm(Species ~ Petal.Width, iris, family = binomial(link = "logit"))
summary(model)
# Valdov test kaze da je Petal.Width znacajan

# Razmotriti kako se tumace koeficijenti

# Kako izgleda nas model...
plot(iris$Petal.Width, as.numeric(iris$Species) - 1)
pl <- seq(1, 2.5, 0.05)
# Prilikom pravljanja predvidjanja, ako nam je potrebna verovatnoca, bitno je
# staviti type="response".
preds <- predict(model, data.frame(Petal.Width = pl), type = "response")
lines(pl, preds)
# Na ovom grafiku su kao tacke prikazani podaci, vrsta moze biti 0 ili 1, dok
# kriva predstavlja ocenu verovatnoce da je vrsta 1 (virginica), u zavisnosti od
# prediktora Petal.Width.

# Posto se bavimo klasifikacijom, hteli bismo da dajemo predvidjanja u obliku 0
# ili 1, umesto verovatnoca. Za to je potrebno samo odrediti prag nakon kojeg
# cemo podatak klasifikovati sa 1. Za pocetak koristicemo prag 0.5, tj. ako je
# verovatnoca da je cvet "virginica" veca od 0.5, klasifikovacemo ga kao
# "virginica", a u suprotnom kao "versicolor".

# Treba nam neki nacin da proverimo koliko dobro klasifikuje model.
# Najjednostavnije je to videti prikazivanjem "matrica konfuzije".
predicted <- predict(model, type="response") > 0.5
confusion <- table(iris$Species, predicted)
# Najbolje je imati sve vrednosti na dijagonali sto znaci da savrseno
# predvidjamo, ali to cesto nije realan zahtev. Takodje moze biti i znak
# preprilagodjavanja (overfitting) trening skupu u kontekstu masinskog ucenja.

# Najpoznatije mere kvaliteta klasifikatora su:
# Oznake: T - "True", F - "False", P - "Positive",  N - "Negative"
# Dakle, TN = True Negative i slicno
# Confusion matrica:
# | TN FP |
# | FN TP |

# Sve ove mere kvaliteta postoje i u paketu Metrics
library(Metrics)
# - Accuracy: udeo tacno predvidjenih podataka
# Accuracy = (TN + TP) / (TN + TP + FN + FP))
sum(diag(confusion)) / sum(confusion) # = 0.94
accuracy(iris$Species == "virginica", predicted) # Metrics paket

# - Precision: udeo pozitivnih podataka od svih proglasenih pozitivnim
# Precision = (TP) / (TP + FP))
confusion[2, 2] / sum(confusion[, 2]) # = 0.96
precision(iris$Species == "virginica", predicted) # Metrics paket

# - Recall: udeo pronadjenih pozitivnih podataka od svih stvarno pozitivnih
# Recall = (TP) / (TP + FN))
confusion[2, 2] / sum(confusion[2, ]) # = 0.92
recall(iris$Species == "virginica", predicted) # Metrics paket

# - F1 score: harmonijska sredina Precision i Recall
prec <- confusion[2, 2] / sum(confusion[, 2])
rec <- confusion[2, 2] / sum(confusion[2, ])
f1 <- 2 * (prec * rec) / (prec + rec) # = 0.94

fbeta_score(iris$Species == "virginica", predicted) # Metrics paket
# Odmera baziranih na matrici konfuzije, f1 score je cesto najbolji za
# koriscenje, jer je otponiji na neizbalansiranost klasa od Accuracy, a uzima u
# obzir i precision i recall.

# Drugi tip mera, koji je otporniji na neizbalansirane klase je ROC kriva i AUC
# - povrsina ispod ROC krive. AUC se moze interpretirati kao da ocenjuje
# verovatnocu da pri slucajnom izboru dve instance iz razlicitih klasa ona iz
# negativne klase ima manji skor od one iz pozitivne klase.

# AUC se dobija funkcijom auc iz Metrics paketa
# NAPOMENA. AUC zahteva da joj se predaju kao drugi argument verovatnoce klasa,
# a ne 0 i 1 vrednosti klasifikacije
Metrics::auc(iris$Species == "virginica", predict(model, type="response"))
# Vrednost od 0.98 je vrlo velika sto znaci da dobro klasifikujemo
# Primetimo da AUC ne zavisi od odabranog praga za klasifikaciju.

# Nas dobijeni model ima prilicno visoke vrednosti svih mera. Najbolju vrednost
# praga mozemo odrediti unakrsnom validacijom ili podelom na trening i
# validacioni skup, sto je jednostavnije i sad cemo uraditi, mada je unakrsna
# validacija preporucljivija.
# Uzecemo stratifikovan uzorak da bismo imali reprezentativniju raspodelu uzorka
set.seed(216)
ind_val1 <- sample(1:50, 20)
ind_val2 <- sample(1:50, 20)
versicolor <- iris[iris$Species == "versicolor", ]
virginica <- iris[iris$Species == "virginica", ]
train <- rbind(versicolor[-ind_val1, ], virginica[-ind_val2, ])
validation <- rbind(versicolor[ind_val1 ,], virginica[ind_val2, ])

m <- glm(Species ~ Petal.Width, train, family = binomial())
preds <- predict(m, validation, type="response")

thresholds <- seq(0.1, 0.9, 0.1)
scores <- sapply(thresholds, function(thr) {
  fbeta_score(validation$Species == "virginica", preds > thr)
})

# najbolji prag...
thresholds[which.max(scores)]
# ispada 0.5
scores
# i za 0.6 i 0.7 je isti f1 skor, pa se mogu uzeti i te vrednosti praga

# Koristeci paket pROC moze se nacrtati ROC kriva, koja opisuje vrednosti "true
# positive rate" (sensitivity) i "false positive rate" (1-specificity) za razne
# vrednosti pragova. Zeljeni oblik grafika je "sto veca grba ka gornjem levom
# uglu", sto indikuje visok nivo pronadjenih pozitivnih instanci i nizak nivo
# instanci pogresno oznacenih kao pozitivne.
library(pROC)
rc <- roc(iris$Species == "virginica", predict(model, type="response"))
plot(rc)
# U nasem slucaju grafik je vrlo povoljan.
# AUC meri povrsinu ispod ove krive.