challenger <- read.delim('challenger.txt') 

logistic_model <- glm(oring ~ temp, data = challenger, family = 'binomial') 
summary(logistic_model)

plot(oring ~ temp, data = challenger)
xn = seq(50, 90, 0.1)
yn = predict(logistic_model, data.frame(temp = xn), type = 'response')
points(xn, yn, type = "l", col = "blue", lwd = 2)

predicted_probabilities = logistic_model$fitted.values
predicted_classes = as.numeric(predicted_probabilities > 0.5)
mean(challenger$oring == predicted_classes) #tacnost modela (accuracy)
confusion_matrix = table(challenger$oring, predicted_classes) #matrica konfuzije
confusion_matrix
# Iz matrice konfuzije možemo da pročitamo razne bitne klasifikacione metrike

# Pomenut Accuracy = (TP + TN)/(P + N) = (TP + TN)/(TP + FP + TN + FN)
acc = (confusion_matrix[1,1] + confusion_matrix[2,2]) / sum(confusion_matrix)
acc

# Sensitivity/Recall/TPR = TP/P = TP/(TP + FN)
sens = confusion_matrix[2,2] / (confusion_matrix[2,1] + confusion_matrix[2,2])
sens

# Specificity/TNR = TN/N = TN/(TN + FP)
spec = confusion_matrix[1,1] / (confusion_matrix[1,1] + confusion_matrix[1,2])
spec

# Precision/PPV = TP/(TP + FP)
prec = confusion_matrix[2,2] / (confusion_matrix[1,2] + confusion_matrix[2,2])
prec

# F1-score: 2*PPV*TPR/(PPV + TPR), harmonijska sredina sens i prec
f1 = 2*prec*sens/(prec + sens)


# Iako je prirodno za prag klasifikacije uzeti C = 0.5, odnosno kada model predviđa verovatnoću p > C klasifikovati pozitivnom klasom,
# u određenim situacijama je poželjno uzeti drugačiji prag za klasifikaciju 
# Na primer, kada je jedna klasa mnogo pristunija od druge, to u modelu pravi prirodnu pristrasnost ka toj klasi
# Što je veći prag C, to model ređe klasifikuje pozitivnom klasom, odnosno teže prepoznaje pozitivnu klasu
# Tada Sensitivity opada, ali Specificity i Precision rastu


# ROC kriva je funkcija Sensitivity od Specificity
library("pROC")
roc_curve = roc(response = challenger$oring, predictor = predicted_probabilities)
plot(roc_curve)
# Želimo da što manje smanjimo Specificity, a postignemo visok Sensitivity
# Zato je model bolji, što je ova kriva viša, odnosno što je površina ispod nje veća
# AUC (Area Under Curve) predstavlja tu metriku
auc = auc(roc_curve)
auc
# Ova površina ima i drugu interpretaciju
# Ona predstavlja statističku meru za to koliko model dobro 'razdvaja' pozitivnu klasu od negativne
# Kao takva ona ne uzima u obzir prag C, i otporna je na nebalansirane podatke


# Logistička regresija sa više atributa
library(mlbench)
data("PimaIndiansDiabetes")
data = PimaIndiansDiabetes
?PimaIndiansDiabetes #informacije o bazi

full_model <- glm(diabetes ~ ., data = data, family = 'binomial')#model sa svim atributima
summary(full_model)

# Predvidjanja radimo standardno:
predicted_probabilities = full_model$fitted.values
predicted_classes = as.numeric(predicted_probabilities > 0.5)

# ROC kriva
roc_curve = roc(response = data$diabetes, predictor = predicted_probabilities)
plot(roc_curve)
auc(roc_curve)

# Matrica konfuzije
confusion_matrix = table(data$diabetes, predicted_classes)
confusion_matrix

# Racunamo tacnost modela
acc = (confusion_matrix[1,1] + confusion_matrix[2,2]) / (confusion_matrix[1,1] + confusion_matrix[1,2] + confusion_matrix[2,1] + confusion_matrix[2,2])
acc #solidna vrednost


# Racunamo odziv modela, posto nam je bitno da za osobe koje stvarno imaju dijabetes, to i predvidimo
sens = confusion_matrix[2,2] / (confusion_matrix[2,1] + confusion_matrix[2,2])
sens #los odziv

# Racunamo i f1-score
prec = confusion_matrix[2,2] / (confusion_matrix[1,2] + confusion_matrix[2,2])
f1 = 2*sens*prec/(sens + prec)
f1

# U ovoj situaciji bismo potencijalno hteli bolji odziv, jer zelimo da za dijabeticare
# model cesce uvidja da imaju dijabetes. Ovo postizemo tako sto promenimo prag klasifikacije.
# Ako zelimo odziv od bar 80%, vidimo sa ROC krive da ce specificnost pasti na oko 70%, sto je u redu.
# Iz podataka za ROC krivu izvlacimo najvece C koje nam to garantuje
new_C = roc_curve$thresholds[max(which(roc_curve$sensitivities > 0.8))]
new_C

new_predicted_classes = as.numeric(predicted_probabilities > new_C)
new_conf_matrix = table(data$diabetes, new_predicted_classes)
new_conf_matrix

acc = (new_conf_matrix[1,1] + new_conf_matrix[2,2]) / (new_conf_matrix[1,1] + new_conf_matrix[1,2] + new_conf_matrix[2,1] + new_conf_matrix[2,2])
acc #tacnost je opala

sens = new_conf_matrix[2,2] / (new_conf_matrix[2,1] + new_conf_matrix[2,2])
sens #80%

prec = new_conf_matrix[2,2] / (new_conf_matrix[1,2] + new_conf_matrix[2,2])
f1 = 2*sens*prec/(sens + prec)
f1 #f1 malo bolji