##### Rozdział 10: Ewaluacja działania modelu -------------------

## Macierze błędów R ----

## Wczytujemy prawdopodobieństwa przewidziane przez klasyfikator SMS-ów z rozdziału 4.
sms_results <- read.csv("sms_results.csv", stringsAsFactors = TRUE)

# kilka pierwszych przypadków testowych
head(sms_results)

# przypadki testowe, w których model jest mniej pewny
head(subset(sms_results, prob_spam > 0.40 & prob_spam < 0.60))

# przypadki testowe, w których model się mylił
head(subset(sms_results, actual_type != predict_type))

# określamy wektory
table(sms_results$actual_type, sms_results$predict_type)

# alternatywne rozwiązanie wykorzystujące interfejs formuł (niepokazane w książce)
xtabs(~ actual_type + predict_type, sms_results)

# używamy funkcji CrossTable
library(gmodels)
CrossTable(sms_results$actual_type, sms_results$predict_type)

# obliczenia dokładności i wskaźnika błędów --
# dokładność
(152 + 1203) / (152 + 1203 + 4 + 31)
# wskaźnik błędów
(4 + 31) / (152 + 1203 + 4 + 31)
# wskaźnik błędów = 1 - dokładność
1 - 0.9748201

## Nie tylko dokładność — inne miary trafności ----
library(caret)
confusionMatrix(sms_results$predict_type, sms_results$actual_type, positive = "spam")

# Statystyka kappa
# przykład z klasyfikatorem SMS-ów
pr_a <- 0.865 + 0.109
pr_a

pr_e <- 0.868 * 0.888 + 0.132 * 0.112
pr_e

k <- (pr_a - pr_e) / (1 - pr_e)
k

# obliczamy statystykę kappa z wykorzystaniem pakietu vcd
library(vcd)
Kappa(table(sms_results$actual_type, sms_results$predict_type))

# obliczamy statystykę kappa z wykorzystaniem pakietu irr
library(irr)
kappa2(sms_results[1:2])

# Współczynnik korelacji Matthewsa (MCC)
# dla klasyfikatora SMS: TN = 1203, FP = 4, FN = 31, TP = 152

# MCC = (TP * TN - FP * FN) /
#         sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN))

# obliczamy ręcznie MCC dla klasyfikatora SMS
(152 * 1203 - 4 * 31) /
  sqrt((152 + 4) * (152 + 31) * (1203 + 4) * (1203 + 31))

# obliczamy MCC za pomocą funkcji mcc() z pakietu mltools
library(mltools)
mcc(sms_results$actual_type, sms_results$predict_type)

# obliczamy MCC z wykorzystaniem korelacji Pearsona
cor(ifelse(sms_results$actual_type == "spam", 1, 0),
    ifelse(sms_results$predict_type == "spam", 1, 0))

# Czułość i swoistość
# przykład z klasyfikatorem SMS-ów
sens <- 152 / (152 + 31)
sens

spec <- 1203 / (1203 + 4)
spec

# przykład wykorzystujący pakiet caret
library(caret)
sensitivity(sms_results$predict_type, sms_results$actual_type, positive = "spam")
specificity(sms_results$predict_type, sms_results$actual_type, negative = "ham")

# Precyzja i przywołanie
prec <- 152 / (152 + 4)
prec

rec <- 152 / (152 + 31)
rec

# przykład wykorzystujący pakiet caret
library(caret)
posPredValue(sms_results$predict_type, sms_results$actual_type, positive = "spam")
sensitivity(sms_results$predict_type, sms_results$actual_type, positive = "spam")

# Miara F
f <- (2 * prec * rec) / (prec + rec)
f

f <- (2 * 152) / (2 * 152 + 4 + 31)
f

## Wizualizacja kompromisów ----
library(pROC)
sms_roc <- roc(sms_results$actual_type, sms_results$prob_spam)

# Krzywa ROC dla naiwnego klasyfikatora bayesowskiego
plot(sms_roc, main = "Krzywa ROC dla filtra spamu SMS ",
     col = "blue", lwd = 2, grid = TRUE, ylab="Czułość", xlab="1–swoistość")


# porównujemy z NN 
sms_results_knn <- read.csv("sms_results_knn.csv")
sms_roc_knn <- roc(sms_results$actual_type, sms_results_knn$p_spam)
plot(sms_roc_knn, col = "red", lwd = 2, add = TRUE)

# obliczamy AUC dla naiwnego klasyfikatora bayesowskiego i kNN
auc(sms_roc)
auc(sms_roc_knn)

## Szacowanie przyszłej trafności ----

# partycjonujemy dane
library(caret)
credit <- read.csv("credit.csv", stringsAsFactors = TRUE)

# Metoda wstrzymywania
# używamy losowych identyfikatorów
random_ids <- order(runif(1000))
credit_train <- credit[random_ids[1:500],]
credit_validate <- credit[random_ids[501:750], ]
credit_test <- credit[random_ids[751:1000], ]

# używamy funkcji z pakietu caret
in_train <- createDataPartition(credit$default, p = 0.75, list = FALSE)
credit_train <- credit[in_train, ]
credit_test <- credit[-in_train, ]

# 10-krotna walidacja krzyżowa
set.seed(123) # gwarantuje powtarzalność wyników
folds <- createFolds(credit$default, k = 10)
str(folds)
credit01_test <- credit[folds$Fold01, ]
credit01_train <- credit[-folds$Fold01, ]

## Automatyzacja 10-krotnej walidacji krzyżowej drzewa decyzyjnego C5.0 za pomocą lapply() ----
library(caret)
library(C50)
library(irr)

credit <- read.csv("credit.csv", stringsAsFactors = TRUE)

set.seed(123)
folds <- createFolds(credit$default, k = 10)

cv_results <- lapply(folds, function(x) {
  credit_train <- credit[-x, ]
  credit_test <- credit[x, ]
  credit_model <- C5.0(default ~ ., data = credit_train)
  credit_pred <- predict(credit_model, credit_test)
  credit_actual <- credit_test$default
  kappa <- kappa2(data.frame(credit_actual, credit_pred))$value
  return(kappa)
})

# badamy wyniki 10 prób
str(cv_results)

# obliczamy średnią statystykę kappa z 10 prób
mean(unlist(cv_results))

# obliczamy odchylenie standardowe z 10 prób
sd(unlist(cv_results))

## Próbkowanie bootstrapowe ----

# obliczamy prawdopodobieństwo, że rekord nie zostanie wybrany
(1 - (1/1000))^1000
(1 - (1/100000))^100000
1 / exp(1) # kiedy n dąży do nieskończoności
