O Problema

# Bibliotecas Necessárias 
library(gdata)
library(dplyr)
library(ggplot2)
library(GGally)
library(caret)
library(C50)
library(e1071)

Impeachment é uma palavra de origem inglesa que significa “impedimento” ou “impugnação”, utilizada como um modelo de processo instaurado contra altas autoridades governamentais acusadas de infringir os seus deveres funcionais. Dizer que ocorreu impeachment ao Presidente da República, significa que este não poderá continuar exercendo funções.

Queremos prever os votos dos deputados referentes ao processo de Impeachment contra a presidência da república, através de dados referentes aos deputados. Iremos utilizar um dataset composto de 540 observações (deputados) e 25 variáveis, onde a variável IMPEACHMENT será nossa variável resposta. São elas:

Os dados para realizar a predição foram extraídos do portal www.qmrepresenta.com.br, onde são elencados votações polêmicas da câmara dos deputados e que repercutiram na mídia nacional. Os temas são:

temas <- read.delim("~/Projetos/DataAnalysis/Classification-Impeachment/dados/temas.csv")
temas$tema
##  [1] Cobrança de cursos em universidades públicas      
##  [2] Tributação serviços de internet                    
##  [3] Terrorismo                                      
##  [4] Infaticídio indígena                              
##  [5] Maioridade 1                                    
##  [6] Maioridade 2                                    
##  [7] Financiamento privado para partidos             
##  [8] Financiamento privado para partidos e candidatos
##  [9] Terceirização                                     
## [10] Distritao                                       
## [11] Transgênico                                      
## [12] Reeleição                                         
## [13] Cota para mulheres legislativo                  
## [14] Pensão                                           
## [15] Seguro Desemprego                               
## [16] Tempo Mandato                                   
## [17] Voto Facultativo                                
## [18] Voto impresso                                   
## [19] Coincidência reeleição                             
## 19 Levels: Cobrança de cursos em universidades públicas ...

Análise Descritiva

# Carregando os dados
deputados_impeachment <- read.csv("~/Projetos/DataAnalysis/Classification-Impeachment/dados/deputados_temas_e_impeachment_v1.1.csv", sep=";")
votacao_impeachment <- read.csv("~/Projetos/DataAnalysis/Classification-Impeachment/dados/imp-votacao.csv")

sum(complete.cases(deputados_impeachment))
## [1] 480

Temos um total de 480 dados completos. O que signicia que 60 não estão completos, ou seja, nem todos os deputados participaram de todas as votações nos temas polêmicos. Mais na frente deveremos analisar se isso será um problema para o nosso modelo.

Vamos olhar variáveis “Tema” para verificar como elas estão distribuidas.

unique(deputados_impeachment$tema_1)
## [1] sim      não       não votou abstenção           art. 17 
## Levels:  abstenção art. 17 não não votou sim

Cada tema pode ser votado da seguinte forma:

unique(deputados_impeachment$IMPEACHMENT)
## [1] SIM   NAO   AUSEN ABST  <NA> 
## Levels: ABST AUSEN NAO SIM

A nossa variável resposta, IMPEACHMENT, pode ser da seguinte forma:

Como a nossa variável resposta é o Impeachment, iremos excluir todas as linhas com NA na variável Impeachment. Como temos apenas 9 deputados votando AUSEN ou ABST também optamos por excluir.

deputados_impeachment <- filter(deputados_impeachment, IMPEACHMENT == "SIM" | IMPEACHMENT == "NAO")
deputados_impeachment <- drop.levels(deputados_impeachment)
d <- deputados_impeachment[, -c(1:5)]

ggpairs(d, upper = "blank", lower = "blank", columns = c(1:4), legends = TRUE) + theme_bw()

ggpairs(d, upper = "blank", lower = "blank", columns = c(5:8), legends = TRUE) + theme_bw()

ggpairs(d, upper = "blank", lower = "blank", columns = c(9:12), legends = TRUE) + theme_bw()

ggpairs(d, upper = "blank", lower = "blank", columns = c(13:16), legends = TRUE) + theme_bw()

ggpairs(d, upper = "blank", lower = "blank", columns = c(17:20), legends = TRUE) + theme_bw()

É possível notar que os temas 12 possui uma distribuição diferente dos demais temas, quase todos o deputados votaram da mesma forma para esse tema, o tema possui pouca variância, portanto será pouco útil para o nosso modelo.

Além disso, os temas 6, 7, 8 e 18 não estão no formato esperado, possuem mais saídas.

Observando mais de perto temos:

unique(deputados_impeachment$tema_6)
## [1] não         sim        não votou  6 não votou            obstção     
## [7] 7 art. 17  não votou   abstção     
## 9 Levels:  6 não votou 7 art. 17 abstção não não votou ... sim
unique(deputados_impeachment$tema_7)
##  [1] sim        não         7 não votou não votou   não votou            
##  [7] obstção      art. 17    7 sim      7 não       abstção     
## 11 Levels:  7 não 7 não votou 7 sim abstção art. 17 ... sim
unique(deputados_impeachment$tema_8)
## [1] sim       não        não votou            não votou abstenção   art. 17  
## Levels:  abstenção art. 17 não não votou não votou sim
unique(deputados_impeachment$tema_18)
## [1] sim       não votou  não        nao votou art. 17   abstenção   im       
## Levels: abstenção art. 17 im não não votou nao votou sim

Acreditamos que essa inconsistência pode prejudicar o modelo final, por essa razão decidimos por excluir os temas 6, 7, 8 e 18.

deputados_impeachment$tema_6 <- NULL
deputados_impeachment$tema_7 <- NULL
deputados_impeachment$tema_8 <- NULL
deputados_impeachment$tema_18 <- NULL

Modelo GBM

Antes de criar qualquer modelo é necessário a criação das partições de treino e teste.

Lembrando que a proporção da variável resposta é:

prop.table(table(deputados_impeachment$IMPEACHMENT))
## 
##       NAO       SIM 
## 0.2727273 0.7272727
set.seed(669)
train_idx = createDataPartition(y = deputados_impeachment$IMPEACHMENT, p=.9, list=FALSE)
train = deputados_impeachment[train_idx,]
test = deputados_impeachment[-train_idx,]
names(train) = names(deputados_impeachment) #adicionando cabeçalho aos dados de treino e test
names(test) = names(deputados_impeachment)

Testando se as proporções foram alteradas de forma significativa

prop.table(table(train$IMPEACHMENT))
## 
##       NAO       SIM 
## 0.2735426 0.7264574
prop.table(table(test$IMPEACHMENT))
## 
##       NAO       SIM 
## 0.2653061 0.7346939

Criando um modelo usando o caret para realizar a validação cruzada, com 10 repetições, utilizando Gradient Boosting Machine (gbm).

Buscando em base de dados externos, iremos incorporar os dados dos deputados que fazem parte da bancada evangélica. Pretendemos dessa forma melhorar o nosso modelo. (Ver https://github.com/nazareno/houseofcunha/blob/master/dados/bancada-evangelica.csv)

bancada_evangelica <- read.table("~/Projetos/DataAnalysis/Classification-Impeachment/dados/bancada_evangelica.csv", header=TRUE, quote="\"")


train$bancada_evangelica <- as.factor(train$id_dep %in% bancada_evangelica$id_dep)
test$bancada_evangelica <- as.factor(test$id_dep %in% bancada_evangelica$id_dep)

Como as variáveis id_dep, nome, deputado são constantes para cada deputado, decidimos por não considerar essas variáveis no modelo. As variáveis que eu considero mais importânte para a criação do primeiro modelo são:

colnames(train[,-c(1,2,3,21)])
##  [1] "partido"            "UF"                 "tema_1"            
##  [4] "tema_2"             "tema_3"             "tema_4"            
##  [7] "tema_5"             "tema_9"             "tema_10"           
## [10] "tema_11"            "tema_12"            "tema_13"           
## [13] "tema_14"            "tema_15"            "tema_16"           
## [16] "tema_17"            "tema_19"            "bancada_evangelica"
ctrl = trainControl(method = "cv", number = 10)
labels = as.factor(train$IMPEACHMENT)
set.seed(669)
gbmFit <- train(x=train[,-c(1,2,3,21)], y=labels,
              method = "gbm",
              verbose = FALSE,
              trControl = ctrl)
gbmFit
## Stochastic Gradient Boosting 
## 
## 446 samples
##  18 predictors
##   2 classes: 'NAO', 'SIM' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## 
## Summary of sample sizes: 401, 401, 402, 401, 401, 401, ... 
## 
## Resampling results across tuning parameters:
## 
##   interaction.depth  n.trees  Accuracy   Kappa      Accuracy SD
##   1                   50      0.8835859  0.6842149  0.04027776 
##   1                  100      0.8856566  0.7029100  0.03246427 
##   1                  150      0.8880303  0.7117129  0.03440202 
##   2                   50      0.8969192  0.7283383  0.02816982 
##   2                  100      0.8923737  0.7250335  0.03288689 
##   2                  150      0.8923737  0.7315816  0.04570975 
##   3                   50      0.8858081  0.7170345  0.04950666 
##   3                  100      0.8923737  0.7313677  0.03776110 
##   3                  150      0.8789394  0.7039250  0.05183548 
##   Kappa SD  
##   0.10578190
##   0.08004063
##   0.08432803
##   0.06852053
##   0.07775726
##   0.10986918
##   0.11802967
##   0.08696600
##   0.11617997
## 
## Tuning parameter 'shrinkage' was held constant at a value of 0.1
## 
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## Accuracy was used to select the optimal model using  the largest value.
## The final values used for the model were n.trees = 50, interaction.depth
##  = 2, shrinkage = 0.1 and n.minobsinnode = 10.

O melhor modelo criado pelo caret para o gbm foi um com o n.trees = 50, interaction.depth = 2, shrinkage = 0.1 and n.minobsinnode = 10. Esse será o modelo que iremos utilizar no teste

test_labels = as.factor(test$IMPEACHMENT)
predictions = predict(gbmFit, newdata= test[,-c(1,2,3,21)])
matrix_confusao <- confusionMatrix(data = predictions, test_labels)

gbm_accuracy <- matrix_confusao$overall[1]
matrix_confusao
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction NAO SIM
##        NAO   9   1
##        SIM   4  35
##                                          
##                Accuracy : 0.898          
##                  95% CI : (0.7777, 0.966)
##     No Information Rate : 0.7347         
##     P-Value [Acc > NIR] : 0.004495       
##                                          
##                   Kappa : 0.7174         
##  Mcnemar's Test P-Value : 0.371093       
##                                          
##             Sensitivity : 0.6923         
##             Specificity : 0.9722         
##          Pos Pred Value : 0.9000         
##          Neg Pred Value : 0.8974         
##              Prevalence : 0.2653         
##          Detection Rate : 0.1837         
##    Detection Prevalence : 0.2041         
##       Balanced Accuracy : 0.8323         
##                                          
##        'Positive' Class : NAO            
## 

Podemos observar, na matriz de confusão, que das 49 instâncias o modelo acertou corretamente 44. Resultando em uma acurácia igual a anterior, 89%.

Modelo Floresta Aleatória

Vamos agora utilizar o caret para encontrar a melhor solução utilizando a Floresta Aleatória como classificador.

# Tranformando para numeric
asNumeric <- function(x) as.numeric(x)
factorsNumeric <- function(d) modifyList(d, lapply(d[, sapply(d, is.factor)],   
                                                   asNumeric))
train2 <- factorsNumeric(train)

ctrl = trainControl(method = "cv", number = 10)
labels = as.factor(train$IMPEACHMENT)
set.seed(669)
rfFit <- train(x=train2[,-c(1,2,3,21)], y=labels,
              method = "C5.0", trControl = ctrl)
rfFit
## C5.0 
## 
## 446 samples
##  18 predictors
##   2 classes: 'NAO', 'SIM' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## 
## Summary of sample sizes: 401, 401, 402, 401, 401, 401, ... 
## 
## Resampling results across tuning parameters:
## 
##   model  winnow  trials  Accuracy   Kappa      Accuracy SD  Kappa SD  
##   rules  FALSE    1      0.8634343  0.6284093  0.03651391   0.09258197
##   rules  FALSE   10      0.8789899  0.6793522  0.04835795   0.12602152
##   rules  FALSE   20      0.8813131  0.6876981  0.04564805   0.11860493
##   rules   TRUE    1      0.8679798  0.6457350  0.04331490   0.11225644
##   rules   TRUE   10      0.8654545  0.6478217  0.04351140   0.10701569
##   rules   TRUE   20      0.8745960  0.6716246  0.03913484   0.09934660
##   tree   FALSE    1      0.8613131  0.6237929  0.05171404   0.13152451
##   tree   FALSE   10      0.8745455  0.6649897  0.04597914   0.13284958
##   tree   FALSE   20      0.8700000  0.6535957  0.03765644   0.09703373
##   tree    TRUE    1      0.8724242  0.6610292  0.03417152   0.08877644
##   tree    TRUE   10      0.8654545  0.6475249  0.03659919   0.08789716
##   tree    TRUE   20      0.8722727  0.6555507  0.03144000   0.08347829
## 
## Accuracy was used to select the optimal model using  the largest value.
## The final values used for the model were trials = 20, model = rules
##  and winnow = FALSE.

O melhor modelo criado pelo caret para o gbm foi um com o trials = 20, model = rules and winnow = FALSE. Esse será o modelo que iremos utilizar no teste

test2 <- factorsNumeric(test)
test_labels = as.factor(test$IMPEACHMENT)
predictions = predict(rfFit, newdata= test2[,-c(1,2,3,21)])
matrix_confusao <- confusionMatrix(data = predictions, test_labels)

rf_accuracy <- matrix_confusao$overall[1]
matrix_confusao
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction NAO SIM
##        NAO  10   1
##        SIM   3  35
##                                          
##                Accuracy : 0.9184         
##                  95% CI : (0.804, 0.9773)
##     No Information Rate : 0.7347         
##     P-Value [Acc > NIR] : 0.001276       
##                                          
##                   Kappa : 0.7798         
##  Mcnemar's Test P-Value : 0.617075       
##                                          
##             Sensitivity : 0.7692         
##             Specificity : 0.9722         
##          Pos Pred Value : 0.9091         
##          Neg Pred Value : 0.9211         
##              Prevalence : 0.2653         
##          Detection Rate : 0.2041         
##    Detection Prevalence : 0.2245         
##       Balanced Accuracy : 0.8707         
##                                          
##        'Positive' Class : NAO            
## 

Podemos observar, na matriz de confusão, que das 49 instâncias o modelo acertou corretamente 45. Resultando em uma acurácia igual a 91%.

Modelo KNN

Vamos agora utilizar o caret para encontrar a melhor solução utilizando o kNN como classificador.

ctrl = trainControl(method = "cv", number = 10)
labels = as.factor(train$IMPEACHMENT)
set.seed(669)
knn_model <- train(train2[,-c(1,2,3,21)], labels,
                 method = "knn",
                 preProcess = c("center","scale"),
                 trControl = ctrl)


knn_model
## k-Nearest Neighbors 
## 
## 446 samples
##  18 predictors
##   2 classes: 'NAO', 'SIM' 
## 
## Pre-processing: centered, scaled 
## Resampling: Cross-Validated (10 fold) 
## 
## Summary of sample sizes: 401, 401, 402, 401, 401, 401, ... 
## 
## Resampling results across tuning parameters:
## 
##   k  Accuracy   Kappa      Accuracy SD  Kappa SD  
##   5  0.8588889  0.6234563  0.04702560   0.12141613
##   7  0.8611111  0.6234034  0.03235215   0.08271836
##   9  0.8655556  0.6330842  0.03604251   0.09807870
## 
## Accuracy was used to select the optimal model using  the largest value.
## The final value used for the model was k = 9.

O melhor modelo criado pelo caret para o KNN foi um com o K = 9. Esse será o modelo que iremos utilizar no teste

test_labels = as.factor(test$IMPEACHMENT)
predictions = predict(knn_model, newdata= test2[,-c(1,2,3,21)])
matrix_confusao <- confusionMatrix(data = predictions, test_labels)

knn_accuracy <- matrix_confusao$overall[1]
matrix_confusao
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction NAO SIM
##        NAO  10   1
##        SIM   3  35
##                                          
##                Accuracy : 0.9184         
##                  95% CI : (0.804, 0.9773)
##     No Information Rate : 0.7347         
##     P-Value [Acc > NIR] : 0.001276       
##                                          
##                   Kappa : 0.7798         
##  Mcnemar's Test P-Value : 0.617075       
##                                          
##             Sensitivity : 0.7692         
##             Specificity : 0.9722         
##          Pos Pred Value : 0.9091         
##          Neg Pred Value : 0.9211         
##              Prevalence : 0.2653         
##          Detection Rate : 0.2041         
##    Detection Prevalence : 0.2245         
##       Balanced Accuracy : 0.8707         
##                                          
##        'Positive' Class : NAO            
## 

Podemos observar, na matriz de confusão, que das 49 instâncias o modelo acertou corretamente 45. Resultando em uma acurácia igual a 91%.

Conclusão

Como para esse problema é tão importante prever os deputados que votam Sim quanto o deputados que votam NÃO, não precisamos ter a preocupação de olhar a Sensitivity ou Specificity. Apenas a acurácia já é o suficiente para esse problema.

Fazendo um comparativo entre os três modelos temos o seguinte gráfico.

toPlot <- data.frame(Modelo = c("KNN", "Floresta", "GBM"), Acuracia = c(knn_accuracy, rf_accuracy, gbm_accuracy))

ggplot(toPlot, aes(x = Modelo, y = Acuracia)) +
  geom_bar(stat="identity") + 
  theme_bw()

Os três modelos apresentaram acurácia bastante similar, sendo o GBM o que teve uma acurácia um pouco menor do que os outros. Para esse problema, o modelo mais adequado é a floresta, pois possui tempo de processamento mais rápido e é tão bom quanto o KNN.