## Statystyka praktyczna w data science (R)
## Rozdział 5. Klasyfikacja
# > (c) 2019 Peter C. Bruce, Andrew Bruce, Peter Gedeck

# Importuje wymagane pakiety R.

library(klaR)
library(MASS)
library(dplyr)
library(ggplot2)
library(FNN)
library(mgcv)
library(rpart)

# Definiuje ścieżki do zestawów danych. Jeżeli przechowujesz dane w innym katalogu, wprowadź poniżej stosowne zmiany.

PSDS_PATH <- file.path(dirname(dirname(getwd())))

loan3000 <- read.csv(file.path(PSDS_PATH, 'dane', 'loan3000.csv'))
loan_data <- read.csv(file.path(PSDS_PATH, 'dane', 'loan_data.csv.gz'))
full_train_set <- read.csv(file.path(PSDS_PATH, 'dane', 'full_train_set.csv.gz'))

# porządkuje zmienną wyniku
loan3000$outcome <- ordered(loan3000$outcome, levels=c('paid off', 'default'))
loan_data$outcome <- ordered(loan_data$outcome, levels=c('paid off', 'default'))
full_train_set$outcome <- ordered(full_train_set$outcome, levels=c('paid off', 'default'))

## Naiwny klasyfikator bayesowski
### Rozwiązanie naiwne

naive_model <- NaiveBayes(outcome ~ purpose_ + home_ + emp_len_, 
                          data = na.omit(loan_data))
naive_model$table

new_loan <- loan_data[147, c('purpose_', 'home_', 'emp_len_')]
row.names(new_loan) <- NULL
new_loan

predict(naive_model, new_loan)

print(predict(naive_model, new_loan))

#### Przykład niedostępny w książce

less_naive <- NaiveBayes(outcome ~ borrower_score + payment_inc_ratio + 
                           purpose_ + home_ + emp_len_, data = loan_data)
less_naive$table[1:2]

stats <- less_naive$table[[1]]
graph <- ggplot(data.frame(borrower_score=c(0,1)), aes(borrower_score)) +
  stat_function(fun = dnorm, color='blue', linetype=1, 
                arg=list(mean=stats[1, 1], sd=stats[1, 2])) +
  stat_function(fun = dnorm, color='red', linetype=2, 
                arg=list(mean=stats[2, 1], sd=stats[2, 2])) +
  labs(y='Prawdopodobieństwo')
graph

## Analiza dyskryminacyjna
### Prosty przykład

loan_lda <- lda(outcome ~ borrower_score + payment_inc_ratio,
                data=loan3000)
loan_lda$scaling

pred <- predict(loan_lda)
print(head(pred$posterior))

#### Rysunek 5.1

# pred <- predict(loan_lda)
# lda_df <- cbind(loan3000, prob_default=pred$posterior[,'default'])

# x <- seq(from=.33, to=.73, length=100)
# y <- seq(from=0, to=20, length=100)
# newdata <- data.frame(borrower_score=x, payment_inc_ratio=y)
# pred <- predict(loan_lda, newdata=newdata)
# lda_df0 <- cbind(newdata, outcome=pred$class)

# graph <- ggplot(data=lda_df, aes(x=borrower_score, y=payment_inc_ratio, color=prob_default)) +
#   geom_point(alpha=.6) +
#   scale_color_gradient2(low='white', high='blue') +
#   scale_x_continuous(expand=c(0,0)) + 
#   scale_y_continuous(expand=c(0,0), lim=c(0, 20)) + 
#   geom_line(data=lda_df0, col='darkgreen', size=2, alpha=.8) +
#   theme_bw()
# graph

# # [Książka]
# graph
# dev.off()

pred <- predict(loan_lda)
lda_df <- cbind(loan3000, prob_default=pred$posterior[,'default'])

center <- 0.5 * (loan_lda$mean[1, ] + loan_lda$mean[2, ])
slope <- -loan_lda$scaling[1] / loan_lda$scaling[2]
intercept = center[2] - center[1] * slope

graph <- ggplot(data=lda_df, aes(x=borrower_score, y=payment_inc_ratio, color=prob_default)) +
  geom_point(alpha=.6) +
  scale_color_gradientn(colors=c('#ca0020', '#f7f7f7', '#0571b0')) +
  scale_x_continuous(expand=c(0,0)) + 
  scale_y_continuous(expand=c(0,0), lim=c(0, 20)) + 
  geom_abline(slope=slope, intercept=intercept, color='darkgreen') +
  theme_bw()

graph

## Regresja logistyczna
### Funkcja odpowiedzi logistycznej i logit

logistic_model <- glm(outcome ~ payment_inc_ratio + purpose_ + 
                        home_ + emp_len_ + borrower_score,
                      data=loan_data, family='binomial')
logistic_model
summary(logistic_model)

p <- seq(from=0.01, to=.99, by=.01)
df <- data.frame(p = p,
                 logit = log(p/(1-p)),
                 odds = p/(1-p))

graph <- ggplot(data=df, aes(x=p, y=logit)) +
  geom_line() +
  labs(x = 'p', y='logit(p)') +
  theme_bw()
graph

### Wartości prognozowane na podstawie regresji logistycznej

pred <- predict(logistic_model)
summary(pred)

prob <- 1/(1 + exp(-pred))
summary(prob)

### Interpretacja współczynników i iloraz szans

graph <- ggplot(data=df, aes(x=logit, y=odds)) +
  geom_line() +
  labs(x='log(ilorazu szans)', y='Iloraz szans') +
  ylim(1, 100) +
  xlim(0, 5) +
  theme_bw()
graph

### Regresja logistyczna z krzywymi sklejanymi

logistic_gam <- gam(outcome ~ s(payment_inc_ratio) + purpose_ + 
                      home_ + emp_len_ + s(borrower_score),
                    data=loan_data, family='binomial')
logistic_gam

### Ocena modelu

terms <- predict(logistic_gam, type='terms')
partial_resid <- resid(logistic_gam) + terms
df <- data.frame(payment_inc_ratio = loan_data[, 'payment_inc_ratio'],
                 terms = terms[, 's(payment_inc_ratio)'],
                 partial_resid = partial_resid[, 's(payment_inc_ratio)'])

graph <- ggplot(df, aes(x=payment_inc_ratio, y=partial_resid, solid = FALSE)) +
  geom_point(shape=46, alpha=0.4) +
  geom_line(aes(x=payment_inc_ratio, y=terms), 
            color='red', alpha=0.5, size=1.5) +
  labs(y='Rezyduum cząstkowe') +
  xlim(0, 25) +
  theme_bw()
graph

df <- data.frame(payment_inc_ratio = loan_data[, 'payment_inc_ratio'],
                 terms = 1/(1 + exp(-terms[, 's(payment_inc_ratio)'])),
                 partial_resid = 1/(1 + exp(-partial_resid[, 's(payment_inc_ratio)'])))

graph <- ggplot(df, aes(x=payment_inc_ratio, y=partial_resid, solid = FALSE)) +
  geom_point(shape=46, alpha=0.4) +
  geom_line(aes(x=payment_inc_ratio, y=terms), 
            color='red', alpha=0.5, size=1.5) +
  labs(y='Rezyduum cząstkowe') +
  xlim(0, 25) +
  theme_bw()
graph

## Ewaluacja modeli klasyfikacji
### Macierz błędów

# Macierz błędów
pred <- predict(logistic_gam, newdata=loan_data)
pred_y <- as.numeric(pred > 0)
true_y <- as.numeric(loan_data$outcome=='default')
true_pos <- (true_y==1) & (pred_y==1)
true_neg <- (true_y==0) & (pred_y==0)
false_pos <- (true_y==0) & (pred_y==1)
false_neg <- (true_y==1) & (pred_y==0)
conf_mat <- matrix(c(sum(true_pos), sum(false_pos),
                     sum(false_neg), sum(true_neg)), 2, 2)
colnames(conf_mat) <- c('Yhat = 1', 'Yhat = 0')
rownames(conf_mat) <- c('Y = 1', 'Y = 0')
conf_mat

### Precyzja, czułość i swoistość

# Precyzja
conf_mat[1, 1] / sum(conf_mat[,1])
# Czułość
conf_mat[1, 1] / sum(conf_mat[1,])
# Swoistość
conf_mat[2, 2] / sum(conf_mat[2,])

### Krzywa ROC

idx <- order(-pred)
recall <- cumsum(true_y[idx] == 1) / sum(true_y == 1)
specificity <- (sum(true_y == 0) - cumsum(true_y[idx] == 0)) / sum(true_y == 0)
roc_df <- data.frame(recall = recall, specificity = specificity)

graph <- ggplot(roc_df, aes(x=specificity, y=recall)) +
  geom_line(color='blue') + 
  scale_x_reverse(expand=c(0, 0)) +
  scale_y_continuous(expand=c(0, 0)) + 
  geom_line(data=data.frame(x=(0:100) / 100), aes(x=x, y=1-x),
            linetype='dotted', color='red') +
  theme_bw() + theme(plot.margin=unit(c(5.5, 10, 5.5, 5.5), "points"))
graph

### Pole pod wykresem krzywej ROC

sum(roc_df$recall[-1] * diff(1-roc_df$specificity))
head(roc_df)

graph <- ggplot(roc_df, aes(specificity)) +
  geom_ribbon(aes(ymin=0, ymax=recall), fill='blue', alpha=.3) +
  scale_x_reverse(expand=c(0, 0)) +
  scale_y_continuous(expand=c(0, 0)) +
  labs(y='Czułość') +
  theme_bw() + theme(plot.margin=unit(c(5.5, 10, 5.5, 5.5), "points"))
graph

## Strategie dla niezbilansowanych danych
### Undersampling

mean(full_train_set$outcome=='default')

full_model <- glm(outcome ~ payment_inc_ratio + purpose_ + home_ + 
                            emp_len_+ dti + revol_bal + revol_util,
                  data=full_train_set, family='binomial')
pred <- predict(full_model)
mean(pred > 0)

 mean(full_train_set$outcome=='default') / mean(pred > 0)

### Oversampling i zwiększenie/obniżenie wag

wt <- ifelse(full_train_set$outcome=='default', 
             1 / mean(full_train_set$outcome == 'default'), 1)
full_model <- glm(outcome ~ payment_inc_ratio + purpose_ + 
                  home_ + emp_len_+ dti + revol_bal + revol_util,
                  data=full_train_set, weight=wt, family='quasibinomial')
pred <- predict(full_model)
mean(pred > 0)

### Generowanie danych¶
#Istnieje wiele różnych implementacji algorytmu SMOTE w R. Pakiet `unbalanced` zawiera algorytm SMOTE i inne metody. Niestety, nie działa on z naszym zestawem danych.
# 
# Na szczęście działa implementacja SMOTE z pakietu `DMwR`.

loan_data_samp <- sample_frac(full_train_set, .05)

# # install.packages('unbalanced')
# library(unbalanced)
# head(full_train_set)
# smote_data <- ubSMOTE(loan_data_samp, loan_data_samp$outcome, 
#                       perc.over = 2000, k = 5, perc.under = 100)
# head(smote_data)

# install.packages('DMwR')
library(DMwR)
smote_data <- SMOTE(outcome ~ ., loan_data_samp, 
                    perc.over = 2000, perc.under=100)
dim(loan_data_samp)
dim(smote_data)
head(smote_data)

### Badanie prognozy

loan_tree <- rpart(outcome ~ borrower_score + payment_inc_ratio,
                   data=loan3000, 
                   control = rpart.control(cp=.005))

geom_abline(slope=slope, intercept=intercept, color='darkgreen')
lda_pred <- data.frame(borrower_score = c((0 - intercept) / slope, (20 - intercept) / slope),
                       payment_inc_ratio = c(0, 20),
                       method = rep('LDA', 2))

tree_pred <- data.frame(borrower_score = c(0.375, 0.375, 0.475, 0.475, 0.575, 0.575),
                        payment_inc_ratio = c(0, 4.426,  4.426, 10.42, 10.42, 20),
                        method = rep('Tree', 6))

glm0 <- glm(outcome ~ (payment_inc_ratio) +  (borrower_score),
            data=loan3000, family='binomial')
y <- seq(from=0, to=20, length=100)
x <- (-glm0$coefficients[1] - glm0$coefficients[2]*y)/glm0$coefficients[3]
glm0_pred <- data.frame(borrower_score=x, payment_inc_ratio=y, method='Logistic')

gam1 <- gam(outcome ~ s(payment_inc_ratio) +  s(borrower_score),
            data=loan3000, family='binomial')
gam_fun <- function(x){
  rss <- sum(predict(gam1, newdata=data.frame(borrower_score=x, payment_inc_ratio=y))^2)
}
est_x <- nlminb(seq(from=.33, to=.73, length=100), gam_fun )
gam1_pred <- data.frame(borrower_score=est_x$par, payment_inc_ratio=y, method="GAM")

loan_fits <- rbind(lda_pred,
                   tree_pred,
                   glm0_pred,
                   gam1_pred)

graph <- ggplot(data=loan_fits, aes(x=borrower_score, y=payment_inc_ratio, color=method, linetype=method)) +
  geom_line(size=1.5) +
  theme(legend.key.width = unit(2,"cm")) +
  guides(linetype = guide_legend(override.aes = list(size = 1)))
graph
