Carregando os pacotes a serem usados.
setwd('/home/heitor/Área de Trabalho/R Projects/Análise Macro/Lab 7')
library(tidyverse) # general
library(rpart) # my standard tree package
library(rpart.plot) # extension to rpart
#library(party) # 1 alternative package
library(tree) # 2 alternative package
library(gt) # to report's tables
library(plotly) # interacting with plotOs dados são informações de antecedentes e características de tomadores de empréstimos e se eles se tornaram ou não inadimplentes. Devemos converter toda a variável em caracteres para fatores, pois elas realmente representam níveis.
dd <- read_csv("credit.csv") %>% as_tibble()
dd[sapply(dd, is.character)] <- lapply(dd[sapply(dd, is.character)], as.factor)
attach(dd)Uma rápida olhada na natureza dos dados:
dd %>% summary()## checking_balance months_loan_duration credit_history
## < 0 DM :274 Min. : 4.0 critical :293
## > 200 DM : 63 1st Qu.:12.0 good :530
## 1 - 200 DM:269 Median :18.0 perfect : 40
## unknown :394 Mean :20.9 poor : 88
## 3rd Qu.:24.0 very good: 49
## Max. :72.0
## purpose amount savings_balance
## business : 97 Min. : 250 < 100 DM :603
## car :337 1st Qu.: 1366 > 1000 DM : 48
## car0 : 12 Median : 2320 100 - 500 DM :103
## education : 59 Mean : 3271 500 - 1000 DM: 63
## furniture/appliances:473 3rd Qu.: 3972 unknown :183
## renovations : 22 Max. :18424
## employment_duration percent_of_income years_at_residence age
## < 1 year :172 Min. :1.000 Min. :1.000 Min. :19.00
## > 7 years :253 1st Qu.:2.000 1st Qu.:2.000 1st Qu.:27.00
## 1 - 4 years:339 Median :3.000 Median :3.000 Median :33.00
## 4 - 7 years:174 Mean :2.973 Mean :2.845 Mean :35.55
## unemployed : 62 3rd Qu.:4.000 3rd Qu.:4.000 3rd Qu.:42.00
## Max. :4.000 Max. :4.000 Max. :75.00
## other_credit housing existing_loans_count job dependents
## bank :139 other:108 Min. :1.000 management:148 Min. :1.000
## none :814 own :713 1st Qu.:1.000 skilled :630 1st Qu.:1.000
## store: 47 rent :179 Median :1.000 unemployed: 22 Median :1.000
## Mean :1.407 unskilled :200 Mean :1.155
## 3rd Qu.:2.000 3rd Qu.:1.000
## Max. :4.000 Max. :2.000
## phone default
## no :596 no :700
## yes:404 yes:300
##
##
##
##
dd %>% glimpse()## Rows: 1,000
## Columns: 17
## $ checking_balance <fct> < 0 DM, 1 - 200 DM, unknown, < 0 DM, < 0 DM, unkn…
## $ months_loan_duration <dbl> 6, 48, 12, 42, 24, 36, 24, 36, 12, 30, 12, 48, 12…
## $ credit_history <fct> critical, good, critical, good, poor, good, good,…
## $ purpose <fct> furniture/appliances, furniture/appliances, educa…
## $ amount <dbl> 1169, 5951, 2096, 7882, 4870, 9055, 2835, 6948, 3…
## $ savings_balance <fct> unknown, < 100 DM, < 100 DM, < 100 DM, < 100 DM, …
## $ employment_duration <fct> > 7 years, 1 - 4 years, 4 - 7 years, 4 - 7 years,…
## $ percent_of_income <dbl> 4, 2, 2, 2, 3, 2, 3, 2, 2, 4, 3, 3, 1, 4, 2, 4, 4…
## $ years_at_residence <dbl> 4, 2, 3, 4, 4, 4, 4, 2, 4, 2, 1, 4, 1, 4, 4, 2, 4…
## $ age <dbl> 67, 22, 49, 45, 53, 35, 53, 35, 61, 28, 25, 24, 2…
## $ other_credit <fct> none, none, none, none, none, none, none, none, n…
## $ housing <fct> own, own, own, other, other, other, own, rent, ow…
## $ existing_loans_count <dbl> 2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2…
## $ job <fct> skilled, skilled, unskilled, skilled, skilled, un…
## $ dependents <dbl> 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ phone <fct> yes, no, no, no, no, yes, no, yes, no, no, no, no…
## $ default <fct> no, yes, no, no, yes, no, no, no, no, yes, yes, y…
gg2 <- ggplot(dd) +
geom_bar(aes(credit_history, fill=default)) +
labs(title = 'Proporção de Defaults por Classificação de Histórico',
x = 'Classificações do Histórico de Empréstimos',
y = 'Contagem')
ggplotly(gg2)gg3 <- ggplot(dd) +
geom_bar(aes(job, fill=default)) +
labs(title = 'Proporção de Defaults por Classificação de Emprego',
x = 'Tipos de Empregos',
y = 'Contagem')
ggplotly(gg3)gg4 <- ggplot(dd) +
geom_bar(aes(purpose, fill=default)) +
labs(title = 'Proporção de Defaults por Aplicação do Empréstimo',
x = 'Destinações',
y = 'Contagem')
ggplotly(gg4)gg5 <- ggplot(dd) +
geom_bar(aes(savings_balance, fill=default)) +
labs(title = 'Proporção de Defaults por Grupos de Poupança',
x = 'Agrupamentos de Poupança',
y = 'Contagem')
ggplotly(gg5)Foi aleatoriamente escolhida uma amostra de treino, sobre o qual o modelo será feito, de 80%. Portanto, o teste será feito na amostra restante, 20%. A aleatoriedade está “fixa” em 777 para fins de replicação.
set.seed(777)
train <- dd %>% slice_sample(.,prop= 0.8)
sid <- as.numeric(rownames(train))
test <- dd[-sid,]
remove(sid)Crontruiremos uma árvore de classificação com os dados de treino, usando todas as demais variáveis como explicativas de default.
tree1 <- rpart(formula = default~.,
data = train,
method = 'class') #"class", classification; "anova", regression
tree1 %>% printcp()##
## Classification tree:
## rpart(formula = default ~ ., data = train, method = "class")
##
## Variables actually used in tree construction:
## [1] amount checking_balance credit_history
## [4] employment_duration months_loan_duration purpose
## [7] savings_balance
##
## Root node error: 234/800 = 0.2925
##
## n= 800
##
## CP nsplit rel error xerror xstd
## 1 0.038462 0 1.00000 1.00000 0.054986
## 2 0.019231 4 0.80342 0.94444 0.054047
## 3 0.015670 6 0.76496 0.96581 0.054419
## 4 0.010000 9 0.71795 0.97863 0.054636
tree1 %>% rpart.plot()tree1 %>% plotcp()Usaremos agora o modelo treinado, tree1 para prever os dados do teste.
pred_tree1 <- predict(tree1,
newdata = test,
type = 'matrix')
pred_tree1 %>% head()## [,1] [,2] [,3] [,4] [,5] [,6]
## 1 1 319 49 0.8668478 0.1331522 0.46000
## 2 1 76 31 0.7102804 0.2897196 0.13375
## 3 1 76 31 0.7102804 0.2897196 0.13375
## 4 1 319 49 0.8668478 0.1331522 0.46000
## 5 1 76 31 0.7102804 0.2897196 0.13375
## 6 2 47 82 0.3643411 0.6356589 0.16125
Transformaremos os resultados previstos em uma Tabela de Confusão, comp_1, de onde conseguimos ver a taxa de acurácia, de falsos negativo e positivo.
comp_1 <- data.frame(Original = test$default,
Previsão = pred_tree1[,1]) %>%
mutate( Previsão = recode( Previsão, '1'='no', '2'='yes'))
comp_1 <- ftable(comp_1)
acc1 <- sum(diag(comp_1))/sum(comp_1)
falsYes1 <- comp_1[1,2]/sum(comp_1[,2])
falsNo1 <- comp_1[2,1]/sum(comp_1[,1])Agora podaremos a primeira árvore, escolhendo uma penalização, CP de acordo com o menor erro relativo \(1-R^2.\), para formar tree2, e com a menor soma dos quadrados dos erros dos resíduos para formar tree3.
tree2 <- prune(tree1,
cp= tree1$cptable[which.min(tree1$cptable[,"rel error"]),"CP"])
tree3 <- prune(tree1,
cp= tree1$cptable[which.min(tree1$cptable[,"xerror"]),"CP"])
# select the complexity parameter associated with the smallest cross-validated errortree2 %>% rpart.plot()tree3 %>% rpart.plot()Usaremos as duas modalidades de poda para fazer a previsão:
pred_tree2 <- predict(tree2,
newdata = test,
type = 'matrix')
pred_tree3 <- predict(tree3,
newdata = test,
type = 'matrix')Faremos agora a Tabela de Confusão para ambas as podas.
comp_2 <- data.frame(Original = test$default,
Previsão = pred_tree2[,1]) %>%
mutate( Previsão = recode( Previsão, '1'='no', '2'='yes'))
comp_3 <- data.frame(Original = test$default,
Previsão = pred_tree3[,1]) %>%
mutate( Previsão = recode( Previsão, '1'='no', '2'='yes'))
comp_2 <- ftable(comp_2)
comp_3 <- ftable(comp_3)
comp_2## Previsão no yes
## Original
## no 123 16
## yes 29 32
comp_3## Previsão no yes
## Original
## no 126 13
## yes 31 30
Arquivaremos os resultados obtidos até agora num data frame novo, compare.
acc2 <- sum(diag(comp_2))/sum(comp_2)
falsYes2 <- comp_2[1,2]/sum(comp_2[,2])
falsNo2 <- comp_2[2,1]/sum(comp_2[,1])
acc3 <- sum(diag(comp_3))/sum(comp_3)
falsYes3 <- comp_3[1,2]/sum(comp_3[,2])
falsNo3 <- comp_3[2,1]/sum(comp_3[,1])
compare <- data.frame(Modelos = c('Tree 1', 'Tree 2', 'Tree 3'),
Accuracy = c(acc1, acc2, acc3),
Falso_No = c(falsNo1, falsNo2, falsNo3),
Falso_Yes = c(falsYes1, falsYes2, falsYes3))
compare %>% gt()| Modelos | Accuracy | Falso_No | Falso_Yes |
|---|---|---|---|
| Tree 1 | 0.775 | 0.1907895 | 0.3333333 |
| Tree 2 | 0.775 | 0.1907895 | 0.3333333 |
| Tree 3 | 0.780 | 0.1974522 | 0.3023256 |
Fiz essa árvore usando outro pacote, para comparar os resultados.
tree4 <- tree(default~.,train)
tree4 %>% summary()##
## Classification tree:
## tree(formula = default ~ ., data = train)
## Variables actually used in tree construction:
## [1] "checking_balance" "amount" "employment_duration"
## [4] "months_loan_duration" "credit_history" "savings_balance"
## [7] "age"
## Number of terminal nodes: 10
## Residual mean deviance: 0.9566 = 755.7 / 790
## Misclassification error rate: 0.235 = 188 / 800
tree4 %>% plot()
tree4 %>% text()cv_t4 <- cv.tree(tree4, FUN = prune.misclass)
cv_t4## $size
## [1] 10 6 5 1
##
## $dev
## [1] 227 227 229 240
##
## $k
## [1] -Inf 0.00 3.00 10.75
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
par(mfrow=c(1,2))
plot(cv_t4$size, cv_t4$dev, type="b")
plot(cv_t4$k, cv_t4$dev, type="b")tree5 <- prune.misclass(tree4, best = cv_t4$size[which.min(cv_t4$dev)])
tree6 <- prune.misclass(tree4, k = cv_t4$k[which.min(cv_t4$dev)])
tree5 %>% summary()##
## Classification tree:
## tree(formula = default ~ ., data = train)
## Variables actually used in tree construction:
## [1] "checking_balance" "amount" "employment_duration"
## [4] "months_loan_duration" "credit_history" "savings_balance"
## [7] "age"
## Number of terminal nodes: 10
## Residual mean deviance: 0.9566 = 755.7 / 790
## Misclassification error rate: 0.235 = 188 / 800
tree6 %>% summary()## Length Class Mode
## size 6 -none- numeric
## dev 1 -none- function
## k 1 -none- numeric
## method 1 -none- character
pred_tree4 <- predict(tree4,test,type="class")
#pred_tree5 <- predict(tree5,test,type="class")
#pred_tree6 <- predict(tree6,test,type="class")comp_4 <- data.frame(Original = test$default,
Previsão = pred_tree4)
comp_4 <- ftable(comp_4)
comp_4## Previsão no yes
## Original
## no 119 20
## yes 27 34
acc4 <- sum(diag(comp_4))/sum(comp_4)
falsYes4 <- comp_4[1,2]/sum(comp_4[,2])
falsNo4 <- comp_4[2,1]/sum(comp_4[,1])
nova_linha <- c('Tree 4', acc4, falsYes4, falsNo4)
compare <- rbind(compare, nova_linha)
remove(acc1, acc2, acc3, acc4,
falsNo1, falsNo2, falsNo3, falsNo4,
falsYes1, falsYes2, falsYes3, falsYes4, nova_linha)
compare %>% gt()| Modelos | Accuracy | Falso_No | Falso_Yes |
|---|---|---|---|
| Tree 1 | 0.775 | 0.190789473684211 | 0.333333333333333 |
| Tree 2 | 0.775 | 0.190789473684211 | 0.333333333333333 |
| Tree 3 | 0.78 | 0.197452229299363 | 0.302325581395349 |
| Tree 4 | 0.765 | 0.37037037037037 | 0.184931506849315 |
Percebemos, portanto, que a melhor árvore é a terceira, pela acurácia geral. Vamos seguir com ela.
Fiz algumas visualizações para tomar mais noção de como estão os dados. Coloquei lado-a-lado Default original e o previsto:
test$Fit_default <- as.factor(pred_tree2[,1])
test <- test %>% mutate( Fit_default = recode( Fit_default, '1'='no', '2'='yes'))
attach(test)## The following objects are masked from dd:
##
## age, amount, checking_balance, credit_history, default, dependents,
## employment_duration, existing_loans_count, housing, job,
## months_loan_duration, other_credit, percent_of_income, phone,
## purpose, savings_balance, years_at_residence
gg1 <- ggplot(test, aes(x=amount, y=months_loan_duration,
colour=default, shape=Fit_default)) +
geom_point() + # poderia colocar o aes em _point, mas, para o lm, precisa ser em gg.
geom_smooth(method=lm, se=FALSE, aes(colour=default,
linetype=Fit_default)) +
labs(title = 'Montante do Empréstimo por Duração',
subtitle = 'Separado por Default e Default Previsto',
x = 'Montante do Empréstimo',
y = 'Duração em Meses')
ggplotly(gg1)## `geom_smooth()` using formula 'y ~ x'