Esta é a minha solução proposta para o desafio da Codenation, para entrar no AceleraDev de Data Science (turma com início no dia 09/06/2020). Mais informações sobre o AceleraDev você encontra aqui:
https://codenation.dev/aceleradev/ds-online-2/
Definição do problema¶
Você deverá criar um modelo para prever a nota da prova de matemática de quem participou do ENEM 2016.
O contexto do desafio gira em torno dos resultados do ENEM 2016 (disponíveis no arquivo train.csv). Este arquivo, e apenas ele, deve ser utilizado para todos os desafios. Qualquer dúvida a respeito das colunas, consulte o Dicionário dos Microdados do Enem 2016.
Muitas universidades brasileiras utilizam o ENEM para selecionar seus futuros alunos e alunas. Isto é feito com uma média ponderada das notas das provas de matemática, ciências da natureza, linguagens e códigos, ciências humanas e redação, com os pesos abaixo:
- matemática: 3
- ciências da natureza: 2
- linguagens e códigos: 1.5
- ciências humanas: 1
- redação: 3
No arquivo test.csv crie um modelo para prever nota da prova de matemática (coluna NU_NOTA_MT) de quem participou do ENEM 2016. Salve sua resposta em um arquivo chamado answer.csv com duas colunas: NU_INSCRICAO e NU_NOTA_MT.
Qualquer dúvida a respeito das colunas, consulte o Dicionário dos Microdados do Enem 2016.
Análise exploratória¶
Pacotes e dados¶
# Pacotes
library(corrplot)
library(dplyr)
library(sqldf)
library(xgboost)
library(methods)
library(ggplot2)
library(Metrics)
library(repr)
library(GGally)
library(gridExtra)
library(MASS)
library(rpart)
library(rpart.plot)
library(randomForest)
library(xgboost)
library(tibble)
library(RColorBrewer)
library(readr)
# Datasets
train = read.csv("./data/train.csv")
test = read.csv("./data/test.csv")
dim(train)
dim(test)
# Nomes das variáveis (colunas) - treino
colnames(train)
colnames(test)
# Pré visualização do conteúdo treino
head(train)
Como temos uma quantidade imensa de dados, fiz a leitura do dicionário de dados, separando os dados em subgrupos para visualizar a estrutura do dataset de treino por completo, como pode ser observado abaixo:
# Dicionário de dados - Treino
data_train_1 = train[, 1:21] # DADOS DO PARTICIPANTE
data_train_2 = train[,22:29] # DADOS DA ESCOLA
data_train_3 = train[,30:42] # DADOS DOS PEDIDOS DE ATENDIMENTO ESPECIALIZADO
data_train_4 = train[,43:47] # DADOS DOS PEDIDOS DE ATENDIMENTO ESPECÍFICO
data_train_5 = train[,48:81] # DADOS DOS PEDIDOS DE RECURSOS ESPECIALIZADOS E ESPECÍFICOS
data_train_6 = train[,82:85] # DADOS DOS PEDIDOS DE CERTIFICAÇÃO DO ENSINO MÉDIO
data_train_7 = train[,86:89] # DADOS DO LOCAL DE APLICAÇÃO DA PROVA
data_train_8 = train[,90:92] # PRESENÇA NAS PROVAS (sem presença na prova de matemática)
data_train_9 = train[93] # PRESENÇA NA PROVA DE MATEMATICA
data_train_10 = train[,94:97] # CÓDIGO DA PROVA
data_train_11 = train[,98:100] # NOTAS DAS PROVAS (sem a nota de matemátia e redação)
data_train_12 = train[,102:110] # VETOR COM GABARITOS
data_train_13 = train[,112:117] # NOTAS DE REDAÇÃO
data_train_14 = train[,118:167] # DADOS DO QUESTIONÁRIO SOCIOECONÔMICO
target_train = train[101] # NOTAS DE MATEMÁTICA
Analisando os datasets, os dados referentes ao desempenho dos alunos encontram-se apenas nas colunas apontadas abaixo. Vou criar um dataset base para a nossa análise do data_train. Como os dois datasets são iguais, com a diferença que o dataset de treino tem menos colunas, farei a análise exploratória apenas com base nos dados de treino.
- Presença na prova de matemática
- Notas das provas
- Notas de redação
# Criando dataset apenas com os dados de desempenho dos alunos
data_train = cbind(data_train_9, data_train_11, data_train_13, target_train)
head(data_train)
Tratamento de dados missing¶
Percebe-se que temos muitos dados nulos, inclusive na nossa variável target. O quanto isso representa no nosso dataset? Vejamos abaixo:
# Função - Proporção dos valores missing
propmiss = function(dataframe) {
m = sapply(dataframe, function(x) {
data.frame(
nmiss = sum(is.na(x)),
n=length(x),
propmiss = sum(is.na(x))/length(x)
)
})
d = data.frame(t(m))
d = sapply(d, unlist)
d = as.data.frame(d)
d$variable = row.names(d)
row.names(d) = NULL
d = cbind(d[ncol(d)],d[-ncol(d)])
return(d[order(d$propmiss), ])
}
# Proporção dos valores missing
propmiss(data_train)
# Proporção de pessoas que faltaram na prova de matemática
(13730 - sum(train$TP_PRESENCA_MT))/13730
Os dados NA aparecem quando o aluno não compareceu no dia da prova. Pela tabela acima, percebemos que 26% dos candidatos faltaram no dia da prova. Como forma de tratamento, podemos substituir os dados missing por 0, considerando que o candidato não tirou nenhuma nota se ele faltou. Precisaremos fazer isso em nosso dataset de teste também, antes de inserir os dados no modelo.
# Removendo valores NA
data_train[is.na(data_train)] = 0
head(data_train)
Correlação dos dados¶
# Plot de correlação
options(repr.plot.width=12, repr.plot.height=9)
m_cor = cor(data_train)
corrplot(m_cor, method = "square", tl.col = "black")
Visualizando o plot de correlação, notamos que as variáveis que mais tem correlação com a nota de matemática são a TP_PRESENCA_MT, referente à presença no dia da prova de matemática, o que é bem óbvio, e a NU_NOTA_LC, referente à prova de linguagens e códigos, que acontece no mesmo dia da prova de matemática.
Já a menor correlação foi com a nota de redação NU_NOTAC_COMP5 mas essa nota apresentou baixa correlação com todas as outras variáveis também, então sem muitas surpresas por aqui.
# Visualizando a maior correlação
options(repr.plot.width=12, repr.plot.height=7)
ggplot(data_train)+
geom_point(aes(x = NU_NOTA_LC, y = NU_NOTA_MT), color = "steelblue")
Parece haver uma relação positiva entre as duas variáveis: quanto maior a NU_NOTA_LC maior a NU_NOTA_MT. Vamos checar o que acontece com a nosso target quando correlacionado com todas as outras variáveis:
# Pairs Graph
options(repr.plot.width=12, repr.plot.height=12)
pairs(data_train[2:11], panel = panel.smooth, col = "steelblue")
O gráfico mostra que existe uma correlação positiva entre todas as variáveis. Quanto maior a nota em uma disciplina, maior a nota em outra disciplina. O que faz muito sentido, imaginando que um aluno que tenha se preparado melhor para matemática, esteja melhor preparado para todas as outras matérias.
Também notamos que há uma concentração de notas zero em praticamente todas as disciplinas, que fogem muito do padrão de concentração das notas médias.
Distribuição das notas¶
# Plot da distribuição das notas
options(repr.plot.width=12, repr.plot.height=7)
ggplot(data = data_train, aes(x = NU_NOTA_MT)) +
geom_histogram(aes(y = ..density..), fill = 'blue', alpha = 0.5,bins = 50) +
geom_density(colour = 'red', lwd = 1) + xlab(expression(bold('Nota obtida'))) +
ylab(expression(bold('Contagem de notas')))
As notas zero, configuram um outlier e isso pode acabar prejudicando a acertividade do nosso modelo. Vejamos o quanto esses número podem influenciar nos resultados:
# Medidas de dispersão (dados completos)
mean(data_train$NU_NOTA_MT) # Média
median(data_train$NU_NOTA_MT) # Mediana
sd(data_train$NU_NOTA_MT) # Desvio padrão
# Medidas de dispersão (Remoção de alunos que faltaram no dia da prova ou que foram eliminados)
data_train2 = data_train[-which(data_train$TP_PRESENCA_MT == 0),] # alunos que faltaram no dia da prova
data_train2 = data_train2[-which(data_train$TP_PRESENCA_MT == 2),] # alunos que foram eliminados da prova
mean(data_train2$NU_NOTA_MT) # Média
median(data_train2$NU_NOTA_MT) # Mediana
sd(data_train2$NU_NOTA_MT) # Desvio padrão
Os valores zerados influenciam em muito nos resultados dos testes, sendo que o desvio padrão foi de 228 para 100. Mas todas as notas zero são referentes à alunos que faltaram ou foram eliminaos? Vamos ver como ficou o histograma agora com as notas removidas.
# Plot da distribuição das notas sem os alunos que faltaram
ggplot(data = data_train2, aes(x = NU_NOTA_MT)) +
geom_histogram(aes(y = ..density..), fill = 'blue', alpha = 0.5,bins = 50) +
geom_density(colour = 'red', lwd = 1) + xlab(expression(bold('Nota obtida'))) +
ylab(expression(bold('Contgem de notas')))
A visualização dos dados ficou muito melhor sem os dias de TP_PRESENCA_MT = 0 e TP_PRESENCA_MT = 2. Temos apenas uma pequena parcela de alunos que compareceram no dia da prova mas que zeraram na disciplina.
Para resolver esse problema (do outlier de zeros), vou seguir com a estratégia de treinar o modelo apenas com os dados de alunos que compareceram. Dados com TP_PRESENCA_MT = 0 e TP_PRESENCA_MT = 2 receberão automaticamente o valor zero.
Preparo dos dados¶
Como falamos anteriormente, utilizarei os dados da presença no dia da prova de matemática para dividir os dados em 2 datasets. Assim, vou usar apenas os dados das notas de alunos que compareceram no dia do teste para criar o dataset de treino.
Dados Teste
# Dados de treino
d1 = train[,98:100] # NOTAS DAS PROVAS (sem a nota de matemátia e redação)
d2 = train[,112:117] # NOTAS DE REDAÇÃO
d3 = train[101] # NOTAS DE MATEMÁTICA
d4 = train[93] # PRESENÇA NA PROVA DE MATEMATICA
# datasets finais para treino
notas = cbind(d1, d2, d3) # Todas as notas
notas[is.na(notas)] = 0 # Substituição de NAs por zero
notas = notas[-which(d4$TP_PRESENCA_MT == 0 | d4$TP_PRESENCA_MT == 2),] # Remoção de faltas e eliminações
Dados validação
Ainda, irei subdividir os resultados em teste e validação, para poder acompanhar as métricas do desempenho dos modelos e escolher o melhor:
# Separação dados de treino e dados de validação
amostra = floor(nrow(notas) * 0.8) # 80% para treino e 20% para teste
notas_train = notas[1:amostra, ] # dataset de treino - primeiros 80%
notas_valid = notas[(amostra + 1):nrow(notas), ] # dataset de teste - 20% restantes
target_valid = notas_valid[[10]] # Vetor com variável target nas notas validação
notas_valid = notas_valid[,1:9] # Dataframe com todas as notas menos matemática
Dados treino
Como não temos a variável TP_PRESENCA_MT no nosso dataset de treino “Linguagens, Códigos e suas Tecnologias", também estará em Matemática, pois são realizados no mesmo dia. Assumirei que todos que não estavam presentes em TP_PRESENCA_LC tiraram zero em Matemática.
# Dados de teste
d1 = test[1] # Identificação do aluno
d2 = test[29:31] # Notas das perguntas multiescolha
d3 = test[34:39] # Notas de redação
data_test = cbind(d1, d2, d3)
data_test[is.na(data_test)] = 0 # Substituição de NAs por zero
# Presença no dia da prova
presenca = test["TP_PRESENCA_LC"] # PRESENÇA NA PROVA DE MATEMATICA = PRESEÇA EM LINGUAGENS E CÓDIGOS
# Separar teste em 2 datasets
data_test0 = data_test[which(presenca$TP_PRESENCA_LC == 0| presenca$TP_PRESENCA_LC == 2),]
data_test1 = data_test[-which(presenca$TP_PRESENCA_LC == 0 | presenca$TP_PRESENCA_LC == 2),]
# Dados para teste
id_test = data_test1[1]
notas_test = data_test1[,2:10]
Uma conferência rápida para ver se os dados de treino teste e validação possuem as mesmas características. Vamos ver se os dados de NU_NOTA_CN (testar com todas) dos três datasets tem a mesma distribuição, ou se certas notas ficaram mais concentradas na hora de separar a validação e o treino:
par(mfcol = c(1,3))
options(repr.plot.width=12, repr.plot.height=5)
hist(notas_train$NU_NOTA_CN, main = "NU_NOTA_CN dos dados de treino", col = "#FB8072")
hist(notas_valid$NU_NOTA_CN, main = "NU_NOTA_CN dos dados de validação", col = "#80B1D3")
hist(notas_test$NU_NOTA_CN, main = "NU_NOTA_CN dos dados de teste", col = "#B3DE69")
OK! Os dados tem mesma distribuição, podemos continuar.
Métricas¶
Para comparar os diversos modelos gerados, resolvi pegar o maior número de métricas possível. Vou utilizar as seguintes:
mae: Mean Absolute Errormdae: Median Absolute Errormse: Mean Squared Errorrae: Relative Absolute Error: sum(ae(actual, predicted)) / sum(ae(actual, mean(actual)))rmse: Root Mean Squared Errorrrse: Root relative squared errorrse: Relative Squared Error : s sse(actual, predicted) / sse(actual, mean(actual))sse: Sum of Squared Errorrsq: r-squared
# Função para o cálculo de r-squared (Não consta no pacote Metrics)
rsq = function(x, y) summary(lm(y~x))$r.squared
# Criando um dataframe vazio para as métricas
metricas = c("nome_incial",0,0,0,0,0,0,0,0,0)
# Função para guardar as metricas de um modelo e retorná-las
calcula_metricas = function(modelo_nome, target_test, ypred) {
# Métricas de desempenho
m1 = round(mae(target_test, ypred),2) # Mean Absolute Error
m2 = round(mdae(target_test, ypred),2) # Median Absolute Error
m3 = round(mse(target_test, ypred),2) # Mean Squared Error
m4 = round(rae(target_test, ypred),2) # Relative Absolute Error: sum(ae(target_test, ypred)) / sum(ae(actual, mean(actual)))
m5 = round(rmse(target_test, ypred),2) # Root Mean Squared Error
m6 = round(rrse(target_test, ypred),2) # root relative squared error
m7 = round(rse(target_test, ypred),2) # Relative Squared Error : s sse(target_test, ypred) / sse(actual, mean(actual))
m8 = round(sse(target_test, ypred),0) # Sum of Squared Error
m9 = round(rsq(target_test, ypred),2) # r-squared
# Vetor de métricas
metricas_do_modelo = c(modelo_nome, m1, m2, m3, m4, m5, m6, m7, m8, m9)
# Inserir as metricas no dataframe de métricas
metricas = rbind(metricas, metricas_do_modelo)
colnames(metricas) = c("nome_modelo", "mae", "mdae", "mse", "rae", "rmse", "rrse", "rse", "sse", "rsq")
metricas <<- metricas
# Salvar as metricas para hisórico
metricas_nome = paste("metricas_", modelo_nome)
assign(metricas_nome, metricas_do_modelo, envir = .GlobalEnv)
}
Modelos¶
Para identificar o melhor modelo para o problema, vou treinar vários para comparar qual deles tem a menor taxa de erros. Começarei pelo mais básico de todos, a regressão linear.
Regressão linear¶
Modelo 1
# Modelo de regressão linear 1 (Com todos os atributos)
modelo_rlin1 = lm(NU_NOTA_MT ~ ., data = notas_train)
summary(modelo_rlin1)
O resumo do modelo acima nos dá a informação de que nem todos os atributos são relevantes para o modelo. Os mais importantes são os marcados com 3 asteriscos NU_NOTA_CN, NU_NOTA_CH, NU_NOTA_LC e NU_NOTA_COMP5e o de importância média NU_NOTA_COMP4. Verificamos assim, que as notas de redação são as que menos influenciam na variável target. Dessa maneira, irei construir um segundo modelo apenas com as variáveis relevante, mas antes, vou salvar os resultados do modelo 1 para fazermos uma comparação.
# Aplicação do modelo nos dados de teste
ypred = predict(modelo_rlin1, notas_valid)
# Métricas de desempenho
calcula_metricas("modelo_rlin1", target_valid, ypred)
metricas = metricas[-1,] # Remover a primeira linha (só da primeira vez)
# Visualização
metricas
Modelo 2
# Modelo de regressão linear 2 (apenas dados relevantes)
modelo_rlin2 = lm(NU_NOTA_MT ~ NU_NOTA_CN + NU_NOTA_CH + NU_NOTA_LC + NU_NOTA_COMP5 , data = notas_train)
summary(modelo_rlin2)
Comparação dos dois modelos
Abaixo, podemos ver que as métricas dos dois modelos são praticamente idênticas. Sendo assim, vou salvar apenas o plot do modelo 2 para fazer a comparação final.
# Aplicação do modelo nos dados de teste
ypred = predict(modelo_rlin2, notas_valid)
# Métricas de desempenho
calcula_metricas("modelo_rlin2", target_valid, ypred)
metricas
# Plot dos resultados
options(repr.plot.width=12, repr.plot.height=7)
predicted = as.data.frame(ypred)
actual = as.data.frame(target_valid)
g1 = ggplot() +
geom_histogram(data = actual, aes(x = target_valid), fill = "darkgray", alpha = 1, bins = 50) +
geom_histogram(data = predicted, aes(x = ypred), fill = "blue", alpha = 0.5, bins = 50) +
ggtitle("Regressão linear 2") + theme(plot.title = element_text(hjust = 0.5))+
xlab("Nota Obtida") + ylab("Contagem das notas")
g1
O Histograma acima mostra como ficaram distribuídas as notas previstas pelo modelo (em azul) em comparação com o dados originais (em cinza). Podemos simplificar afirmando que quanto mais parecido o histograma azul for do cinza, maior a acurácia? O histograma em conjunto com as métricas nos dará uma boa idéia da precisão do modelo.
Árvore de decisão¶
Modelo 1
A árvore de decisão geralmente não é indicada para a criação de modelos com target contínuo. Vamos ver como ela se sai em comparação com os outros modelos. Primeiro, vou fazer um modelo simples com os parâmetros default:
# Árvore de decisão 1
modelo_tree1 = rpart(formula = (NU_NOTA_MT) ~ ., data = notas_train, method = "anova")
# Aplicação do modelo nos dados de teste
ypred = predict(modelo_tree1, notas_valid)
# Métricas de desempenho
calcula_metricas("modelo_tree1", target_valid, ypred)
metricas
Pelas métricas, o modelo ficou muito parecido com o modelo de regressão linear. Vamos ver como o plot ficou com esse modelo:
# Plot dos resultados
options(repr.plot.width=12, repr.plot.height=7)
predicted = as.data.frame(ypred)
actual = as.data.frame(target_valid)
g2 = ggplot() +
geom_histogram(data = actual, aes(x = target_valid), fill = "darkgray", alpha = 01, bins = 50) +
geom_histogram(data = predicted, aes(x = ypred), fill = "blue", alpha = 0.5, bins = 50) +
ggtitle("Árvore de decisão 1") + theme(plot.title = element_text(hjust = 0.5))+
xlab("Nota Obtida") + ylab("Contagem das notas")
g2
Esse modelo é péssimo!! Ele pode apresentar métricas OK, porém, podemos ver que os resultados estão concentrados em pontos específicos. Quais as decisões que construíram essa árvore?
# Plot da árvore
rpart.plot(modelo_tree1)
Aqui está o motivo da decision tree não ser a melhor escolha para variáveis contínuas. Ela possui um número limitado de respostas e pudemos ver graficamente como isso ocorre. Ainda assim, podemos melhorar um pouco esse modelo. O que acontece se mudamos um alguns parâmetros?
Modelo 2
# Parâmetros
control = rpart.control(minbucket = 50, cp = 0.00001)
minbucket: nº míminmo de leaf nodes. No outro modelo tivemos poucas leaf nodes. Vamos aumentar de 7 (defeault) para 50, que é o número de bins que estamos mostrando no histograma.cp: complexity parameter. Default = 0.01 significa que se no próximo split não melhorar o nosso r-squared em pelo menos 0.01, então o algorítimo deve parar. Podemos colocar um valor muito baixo de cp, assim, o algorítimo irá demorar muito para alcançar esse valor e fará mais splits, o que melhorará a cara do nosso gráfico. Vou mudar de 0.01 para 0.0001
# Criação do modelo (com parâmetros de controle)
modelo_tree2 = rpart(formula = (NU_NOTA_MT) ~ ., data = notas_train, method = "anova", control = control)
# Aplicação do modelo nos dados de teste
ypred = predict(modelo_tree2, notas_valid)
# Métricas de desempenho
calcula_metricas("modelo_tree2", target_valid, ypred)
metricas
# Plot dos resultados
options(repr.plot.width=12, repr.plot.height=7)
predicted = as.data.frame(ypred)
actual = as.data.frame(target_valid)
g3 = ggplot() +
geom_histogram(data = actual, aes(x = target_valid), fill = "darkgray", alpha = 01, bins = 50) +
geom_histogram(data = predicted, aes(x = ypred), fill = "blue", alpha = 0.5, bins = 50) +
ggtitle("Árvore de decisão 2") + theme(plot.title = element_text(hjust = 0.5))+
xlab("Nota Obtida") + ylab("Contagem das notas")
g3
As métricas melhoraram um pouco e o gráfico agora parece algo decente. A decision tree pode até não ser a mais indicada quando se trata de dados contínuos mas nesse caso ela se saiu melhor do que a simples regressão linear.
Plot da decision tree
# Plot sem prune
rpart.plot(modelo_tree2)
Resultados melhores... Mas uma árvore monstruosa! Vamos seguir em frente com outros tipos de modelos.
Random Forest¶
O Random Forest é um algorítimo ótimo, porém, o modelo geralmente demora para ser treinado. Vou criar um modelo com os parâmetros default e ver o que acontece:
# Criação do modelo
modelo_rf = randomForest(NU_NOTA_MT ~ ., data=notas_train, proximity=TRUE)
# Resumo do modelo
modelo_rf
# Aplicação do modelo nos dados de teste
ypred = predict(modelo_rf, notas_valid)
# Métricas de desempenho
calcula_metricas("modelo_rf", target_valid, ypred)
metricas
# Plot dos resultados
options(repr.plot.width=12, repr.plot.height=7)
predicted = as.data.frame(ypred)
actual = as.data.frame(target_valid)
g4 = ggplot() +
geom_histogram(data = actual, aes(x = target_valid), fill = "darkgray", alpha = 01, bins = 50) +
geom_histogram(data = predicted, aes(x = ypred), fill = "blue", alpha = 0.5, bins = 50) +
ggtitle("Random Forest") + theme(plot.title = element_text(hjust = 0.5))+
xlab("Nota Obtida") + ylab("Contagem das notas")
g4
As métricas melhoraram um pouco mas o modelo passou de menos de 1 segundo para 1 minuto e meio para ser treinado.
Fiz um teste aumentando o ntrees da Random Forest e o resultado não melhorou, portanto, manteremos apenas o modelo com parâmetros default para a comparação final
O próximo modelo é mais eficiente em termos de velocidade e costuma apresentar uma boa precisão também.
XGBoost¶
Modelo 1
O XGBoost é um algorítimo rápido e que tem uma precisão boa. Vamos treinar com os parâmetros default:
# Separando dados para treino e dados da variável target
data = as.matrix(notas_train[1:9])
label = notas_train[[10]]
# Criando o objeto de entrada - Classe 'xgb.DMatrix'
xgmat = xgb.DMatrix(data, label = label)
str(xgmat)
# Criação do modelo
modelo_xgb1 = xgboost(data = xgmat, nround = 150, verbose = FALSE)
# Dados de entrada (Treino)
data = as.matrix(notas_valid)
xgmat = xgb.DMatrix(data)
# Aplicação do modelo nos dados de teste
ypred = predict(modelo_xgb1, xgmat)
# Métricas de desempenho
calcula_metricas("modelo_xgb1", target_valid, ypred)
metricas
# Plot dos resultados
options(repr.plot.width=12, repr.plot.height=7)
predicted = as.data.frame(ypred)
actual = as.data.frame(target_valid)
g5 = ggplot() +
geom_histogram(data = actual, aes(x = target_valid), fill = "darkgray", alpha = 01, bins = 50) +
geom_histogram(data = predicted, aes(x = ypred), fill = "blue", alpha = 0.5, bins = 50) +
ggtitle("XGBoost 1") + theme(plot.title = element_text(hjust = 0.5))+
xlab("Nota Obtida") + ylab("Contagem das notas")
g5
Em relação às métricas. O XGBoost demonstrou uma precisão menor do que o Random Forest e a Decision Tree 2, e o histograma também ficou bem parecido. Entretando, vale ressaltar que o XGBoost foi extremamente rápido e demorou menos de 1 segundo para treinar o algorítimo, sendo que o Random Forest demorou cerca de 1,5 minutos. Vou tentar dar uma melhorada nesse modelo, como ele é rápido para ser treinado, não custa fazer alguns testes.
Modelo 2
# Parâmetros
param = list("objective" = "reg:gamma")
# Separando dados para treino e dados da variável target
data = as.matrix(notas_train[1:9])
label = notas_train[[10]]
xgmat = xgb.DMatrix(data, label = label)
modelo_xgb2 = xgboost(data = xgmat, nround = 37, params = param ,verbose = FALSE)
data = as.matrix(notas_valid)
xgmat = xgb.DMatrix(data)
ypred = predict(modelo_xgb2, xgmat)
calcula_metricas("modelo_xgb2", target_valid, ypred)
metricas
# Plot dos resultados
options(repr.plot.width=12, repr.plot.height=7)
predicted = as.data.frame(ypred)
actual = as.data.frame(target_valid)
g6 = ggplot() +
geom_histogram(data = actual, aes(x = target_valid), fill = "darkgray", alpha = 01, bins = 50) +
geom_histogram(data = predicted, aes(x = ypred), fill = "blue", alpha = 0.5, bins = 50) +
ggtitle("XGBoost 2") + theme(plot.title = element_text(hjust = 0.5))+
xlab("Nota Obtida") + ylab("Contagem das notas")
g6
Depois de fazer vários testes (que foram omitidos aqui), o melhor modelo que consegui extrair foi esse mostrado acima, onde eu mudei o objective para reg:gamma e diminui o nround de 150 para 37. Com isso, o r-squared aumentou, ficando acima inclusive do Random Forest (por pouco!) e os erros também diminuiram.
Comparação dos modelos¶
Histogramas
Pelos modelos acima, percebemos que o random forest foi melhor em termos de precisão inicialmente, mas o 2º modelo de XGBoost, com parâmetros modificados acabou tendo uma performance um pouco superior. Como isso se reflete graficamente?
Abaixo podemos ver como os histogramas dos modelos (azul) ficaram em comparação com os dados target (cinza) e em comparação com os outros modelos. Podemos verificar se a nossa hipótese de que "quanto mais parecido o histograma azul for do cinza, maior a acurácia" é verdadeira:
grid.arrange(g1 + xlab(NULL) + ylab(NULL),
g2 + xlab(NULL) + ylab(NULL),
g3 + xlab(NULL) + ylab(NULL),
g4 + xlab(NULL) + ylab(NULL),
g5 + xlab(NULL) + ylab(NULL),
g6 + xlab(NULL) + ylab(NULL),
ncol = 3)
Percebe-se que a forma do histograma, nesse caso, onde a precisão do modelo foi relativamente baixa, não reflete o resultado real. Se eu tivesse que definir o melhor modelo com base na forma do histograma, diria que a Regressão linear 2 ou o XGBoost 1 foram os mais precisos pois a distribuição da frequência das notas predicted está mais próxima da dos resultados target. Entretanto, é muito possível que a distribuição seja parecida mas que o valor predicted para cada elemento actual tenha sido discrepante (parece que foi isso mesmo que ocorreu) . Vejamos uma comparação das métricas de desempenho.
Métricas de desempenho
Temos o dataframe de métricas abaixo, que nos confirma o melhor desempenho dos resultados do XGBoost 2. Porém, além de serem muitas métricas (não havia a necessidade de tantas, eu sei), não é possível identificar o melhor modelo de imediato, podemos fazer isso graficamente também. Para isso, vou reduzir a quantidade de métricas pegando apenas as mais usuárias: mae, rmse e rsq.
# Passos para transformação de matriz em dataframe
options(digits=10)
m1 = as.data.frame(metricas, stringsAsFactors = FALSE)
m2 = as.data.frame(sapply(m1[,2:10], as.numeric))
m3 = add_column(m2, m1[,1], .before = "mae")
colnames(m3)[1] = "nome_modelo"
rownames(m3) = NULL
metricas_df = m3
metricas_df
options(repr.plot.width=12, repr.plot.height=9)
g6 = ggplot(metricas_df, aes(x = nome_modelo)) +
geom_col(aes(y = mae, fill = nome_modelo)) + scale_fill_brewer(palette="Set3") + xlab(NULL) + ylab(NULL) +
theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) + ggtitle("Mean Absolute Error (MAE)")+
coord_cartesian(ylim = c(57, 63))
g7 = ggplot(metricas_df, aes(x = nome_modelo)) +
geom_col(aes(y = rmse, fill = nome_modelo)) + scale_fill_brewer(palette="Set3") + xlab(NULL) + ylab(NULL) +
theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) + ggtitle("Root Mean Squared Error (RMSE) ")+
coord_cartesian(ylim = c(74, 80))
g8 = ggplot(metricas_df, aes(x = nome_modelo)) +
geom_col(aes(y = rsq, fill = nome_modelo)) + scale_fill_brewer(palette="Set3") + xlab(NULL) + ylab(NULL) +
theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) + ggtitle("R-Squared (RSQ)")+
coord_cartesian(ylim = c(0.35, 0.45))
grid.arrange(g6, g7, g8, nrow = 3 )
Não se engane, os gráficos acima tiveram seus limites alterados pois os valores eram muito parecidos e não dava para enxergar a diferença direito! Sabendo disso, ainda assim podemos ver que o XGBoost 2 teve um desempenho melhor em comparação com os outros modelos.
O mae é a média de todos os erros absolutos. Ou seja, é a soma da diferença de cada valor actual menos o seu respectivo valor predicted (valores absolutos), dividido pelo número total de elementos. Sendo assim, é uma forma de se medir o erro e quanto menor o erro, melhor o modelo. Pelo gráfico percebemos que o XGBoost 2 teve o menor erro, seguido pelo Random Forest. Já o modelo com pior desempenho nesse quesito foi o XGBoost 1.
O rmse é a raiz quadrada da média dos erros ao quadrado. Também é uma medida de erro e segue um padrão semelhante ao mae.
O rsq, também chamado de coeficiente de determinação, é uma medida estatística de quão próximos os dados estão da linha de regressão ajustada, é a porcentagem da variação da variável resposta que é explicada por um modelo e vai de 0 à 1. Sendo assim, é uma medida de acurácia, e quanto maior o seu valor, melhor o modelo. No caso, novamente, O modelo do XGBoost mostrou-se melhor, seguido pela Random Forest, mas agora, os piores modelos foram os de Regressão Linear..
Conclusão¶
O modelo será julgado pela sua precisão e portanto, utilizarei o XGBoost 2 para enviar os resultados. Vale lembrar que não foi feita nenhuma transformação nos dados: normalização ou divisão dos alunos em clusters. Isso poderia melhorar a precisão do modelo e pode ser uma matéria de estudo de projetos futuros.
Enfim, vamos enviar o arquivo com a resposta!
Salvando resposta
# Separando dados para treino e dados da variável target
data = as.matrix(notas_test)
xgmat = xgb.DMatrix(data)
ypred = predict(modelo_xgb2, xgmat)
# Chamando o modelo Random Forest para os dados de teste
ypred = predict(modelo_xgb2, xgmat)
predicted = as.data.frame(ypred)
# Separar teste em 2 datasets
resposta_0 = data_test0[1] # Numero de inscrição dos alunos faltantes e eliminados
resposta_0$NU_NOTA_MT = 0
resposta_1 = id_test
resposta_1$NU_NOTA_MT = predicted$ypred
# Resposta final
options(digits = 5)
resposta = rbind(resposta_0, resposta_1)
resposta = resposta[order(as.numeric(rownames(resposta))),,drop=FALSE]
write_csv(resposta, "./data/answer.csv")
