# Построение дерева решений для задачи классификации на основе предложенного набора данных. 
# Оценка качества классификации, формирование прогнозной модели в R 
# Работа с пакетами rpart,  caret

library(caTools)     # работа с данными (в том числе, разделение на тестовую и обучающую выборки)
library(rpart)
library(rpart.plot)
library(dplyr)       #mutate
library(caret)       #confusionMatrix
library(MLmetrics)   # расчет метрик
library(DMwR2)
library(randomForest)

options(scipen = 999)

data <-  read.csv("diabetes.csv")
summary(data)

data <-  mutate(data, Outcome = as.factor(Outcome))

set.seed(1) 
split <- sample.split(data, SplitRatio = 0.7)
# SplitRatio - в каком соотношении будет поделена выборка на train/test 

train <- subset(data, split == TRUE)
test  <- subset(data, split == FALSE)

tree_res <-  rpart(Outcome ~ ., data=train, method="class", 
                   control = rpart.control(minsplot=10, minbucket=5, maxdepth=6))
# control правила останова
# maxdepth количество слоев дерева
# minbucket минимальный размер потомка сколько наблюдений может быть минимально в узле потомке
# minsplot минимальное количество в узле родителя, например, minsplot=9 дальше узел не делим

plot(tree_res)
text(tree_res, use.n=T)
summary(train)


rpart.plot(tree_res, type=1)
rpart.plot(tree_res, type=2, extra = 109)
rpart.plot(tree_res, type=3)
rpart.plot(tree_res, type=4,extra = 109, fallen.leaves=F)

print(tree_res, digits=2)

# node), split, n, loss, yval, (yprob)
# node - номер узла. Номера не идут подряд, а следуют правилу 
# "узлы-потомки узла номер n имеют номера 2*n и 2*n+1"
# split - условие, которое должно быть выполнено,
# чтобы наблюдение попало в узел при расщепления узла-родителя
# n - количество наблюдений из обучающей выборки, попавших в узел
# loss - число наблюдений другого класса, попавших в узел (число ошибок)
# yval - какой класс приписывается наблюдению из данного узла. Определяется тем,
# наблюдения какого класса составляют в узле большинство.
# (yprob) - доля объектов из каждого класса в узле. Интерпретируется как вероятность того,
# что объект, попавший в узел, будет принадлежать классу.


predict_train <- predict(tree_res, train)
predict_train <- predict(tree_res, train, type="class")
table(train$Outcome, predict_train)
#   predict_train
#        0   1
#    0 304  27   факт
#    1  66 116   факт

table(train$Outcome)
table(predict_train)


predict_test <- predict(tree_res, test, type="class")
table(test$Outcome, predict_test)
table(test$Outcome)
table(predict_test)

# Матрица смежности (confusion matrix)
matrix_S <- table(test$Outcome, predict_test)
# Доля верно классифицированных объектов
accuracy <- sum(diag(matrix_S)) / sum(matrix_S)
accuracy




#Матрица сопряженности:
#  - data - прогнозные значения
#  - reference - фактические значения
#  - positive - указать какой уровень считаем положительным
#  - mode - указываем "everything", чтобы получить все метрики

confusionMatrix(data      = predict_test, 
                reference = test$Outcome, 
                positive  = "0", mode = "everything")

confusionMatrix(data      = predict_test, 
                reference = test$Outcome, 
                positive  = "1", mode = "everything")

# Accuracy (Точность) - доля всех правильных предсказаний. 
# Mcnemar's Test P-Value - p-значение для статистического теста, сравнивающего количество
# ложноположительных и ложноотрицательных результатов. p-value>0,05 указывает на то, что нет 
# существенной разницы между количеством ложноположительных и ложноотрицательных результатов 
# Sensitivity (Чувствительность) / Recall (Полнота) - доля истинно положительных предсказаний 
# среди всех положительных действительно существующих экземпляров 
# Specificity (Специфичность) - доля истинно отрицательных предсказаний среди 
# всех отрицательных действительно существующих экземпляров  
# Precision - доля истинно положительных предсказаний среди всех предсказанных положительных экземпляров 194/(194+43)
# F1 Score- гармоническое среднее между точностью и полнотой, которое учитывает как ложные положительные, 
# так и ложные отрицательные предсказания 
# Balanced Accuracy: среднее значение чувствительности и специфичности, 
# что обеспечивает сбалансированную оценку эффективности модели в обоих классах.

predict_test <- predict(tree_res, test, type="class")
matrix_S <- table(test$Outcome, predict_test)
matrix_S 
# Доля верно классифицированных объектов
accuracy <- sum(diag(matrix_S)) / sum(matrix_S)
accuracy                    
# Кросс-валидация (K-folds Cross-Validation)
# Метод кросс-валидации позволяет получить оптимальные параметры для построения дерева.
# Разобьем тренировочные данные на k частей (folds). Используем k-1 часть для тренировки модели
# со всеми возможными значениями подбираемого параметра и протестируем получившиеся модели на оставшейся части.

# Имея тестовые данные определим значение параметра, дающее наилучшие результаты в среднем случае. 

fitControl <- trainControl(method="cv", number=30)
# method="cv" - использовать кросс-валидацию
# number=30 - использовать 30 частей(folds)


# Использовть значения cp от 0.001 до 0.2
cartGrid <- expand.grid(.cp=(1:200)*0.001)


# !ВАЖНО: Зависимая переменная должна иметь тип factor
# Кросс-валидация
trainControl <- train(Outcome ~ ., data=train,
                      method="rpart", trControl=fitControl,tuneGrid=cartGrid)   ###Долная процедура
trainControl$results
# Результатом будет таблица, содержащая различные значения точности(Accuracy) для различных cp
# Нам нужен cp, который дает максимальную точность, он указан в конце вывода функции

best <- trainControl$results[which.max(trainControl$results[,"Accuracy"]),"cp"]
best
# Полученный cp необходимо использовать в функции rpart в качестве параметра
# control=rpart.control(cp=)  

tree_res2 <- rpart(Outcome ~ ., data=train, method="class",
                      control=rpart.control(cp = best))

# Матрица смежности (confusion matrix)
predict_test_2 <- predict(tree_res2, test, type="class")
matrix_S_2 <- table(test$Outcome, predict_test_2)
matrix_S_2
# Доля верно классифицированных объектов
accuracy_1 <- sum(diag(matrix_S_2)) / sum(matrix_S_2)
accuracy_1
rpart.plot(credit_res2, type=4,extra = 109, fallen.leaves=F)

# Оценка важности признаков
importance <- tree_res2$variable.importance
importance
# Визуализация важности признаков
barplot(importance, main = "Важность признаков", 
        col = "cornsilk4", las = 2, cex.names = 0.8)

### Модель случайного леса (RandomForest)

# ntree - число деревьев (рекомендация состоит в том, чтобы взять избыточное количество деревьев
# а потом мы с той информацией,  которую выдаст на процедура в результате своей работы
# понять сколько же на самом деле нужно деревьев оставить
# то есть сначала избыточное количество а потом разумным образом
# количество деревьев сократить) 

# mtry число переменных подвыборки (= корню из количества регрессоров)

# sampsize число наблюдений в подвыборке (у Бреймана (1-1/exp)=0.632, в некоторых других источниках 2/3) 

# nodesize минимальное число наблюдений в узле деревьев (от 1 до 10) 

# replace можно ли выбирать одну подвыборку одни и те же строчки, т.е. могут ли в подвыборке быть дубли
# утверждается, что когда можно выбирать дубли как будто результат получше
set.seed(111)
n.tree = 500
rf_res_0 <-  randomForest(train$Outcome ~., data=train, do.trace = n.tree/10)
mtry=floor(sqrt(ncol(train)-1))
set.seed(222)

nodesize_1 = 1
keep.forest_1 = TRUE #записать построенные случайный лес
mtry=4
# do.trace = n.tree/10 показывать промежуточный результат каждые 10% "работы"
rf_res <-  randomForest(Outcome ~ ., data=train, ntree=n.tree, mtry=floor(sqrt(ncol(train)-1)), 
                        replace=TRUE, nodesize = nodesize_1, 
                        importance=TRUE, localImp=FALSE,
                        proxomoty=FALSE, norm.votes = TRUE, 
                        do.trace = n.tree/10, keep.forest = keep.forest_1, 
                        corr.bias = FALSE, keep.inbag=FALSE)
# ntree      OOB      1      2
# 50:      26.71% 19.03% 40.66%
# 100:     27.10% 18.13% 43.41%
# 150:     26.32% 17.52% 42.31%
# 200:     25.93% 16.62% 42.86%
# 250:     24.76% 16.01% 40.66%
# 300:     25.73% 16.62% 42.31%
# 350:     25.15% 16.62% 40.66%
# 400:     24.17% 15.41% 40.11%
# 450:     24.37% 15.41% 40.66%
# 500:     25.15% 16.31% 41.21%
# процент ошибок при распозновании "1 класса" и "2 класса"
# OOB out of back процент ошибок на часть выборки, которая не попала в обучение


cartGrid <- expand.grid(mtry = seq(2, 16, by = 1)) # определение набора параметров для перебора


trainControl <- trainControl(method = "cv", number = 10) # подготовка тренировочной выборки

# Подбор параметров случайного леса
set.seed(222)
rf_tune <- train(Outcome ~ ., data=train, method = "rf", tuneGrid = cartGrid, ntree = 500, trControl = trainControl) #!!!ДОЛГАЯ ПРОЦЕДУРА
print(rf_tune)

# mtry  Accuracy   Kappa    
# 2    0.7605339  0.4625149
# 3    0.7509156  0.4435381
# 4    0.7410377  0.4234856
# 5    0.7468461  0.4355975
# 6    0.7507677  0.4465487
# 7    0.7429985  0.4272769
# 8    0.7487692  0.4412216
# 9    0.7487329  0.4394489
# 10    0.7449593  0.4317929
# 11    0.7507300  0.4435470
# 12    0.7449230  0.4283999
# 13    0.7429622  0.4256993
# 14    0.7390770  0.4211208
# 15    0.7391147  0.4156928
# 16    0.7391147  0.4176845
# 
# Accuracy was used to select the optimal model using the largest value.
# The final value used for the model was mtry = 2.

rf_res <-  randomForest(Outcome~., data=train, ntree=500, mtry=2, 
                          replace=TRUE, nodesize = nodesize_1, 
                          importance=TRUE, localImp=FALSE,
                          proxomoty=FALSE, norm.votes = TRUE, 
                          do.trace = n.tree/10, keep.forest = keep.forest_1, 
                          corr.bias = FALSE, keep.inbag=FALSE)

imp <- importance(rf_res, type=1)
featureImportance <- data.frame(Feature=row.names(imp), Importance=imp[,1])

ggplot(featureImportance, aes(x=reorder(Feature, Importance), y=Importance)) +
  geom_bar(stat="identity", fill="cornsilk4") +
  coord_flip() + 
  theme_light(base_size=16) +
  xlab("Importance") +
  ylab("") + 
  ggtitle("Random Forest Feature Importance\n") +
  theme(plot.title=element_text(size=14))


Y_pred_test <- predict(rf_res, newdata = test, n.trees = 2500)

confusionMatrix(data      = Y_pred_test, 
                reference = test$Outcome, 
                positive  = "0", mode = "everything")