Aprendizado Estatístico de Máquina
Rafael Izbicki
Material explicativo: Machine Learning sob a ótica estatística: Uma abordagem preditivista para estatística com exemplos em R (http://www.rizbicki.ufscar.br/sml)
set.seed(1)
library(data.table)
library(tm)
library(glmnet)
library(tidyverse)
dados <- fread("../../../../mineracaoDeDados/2016_01/notasDeAula/data/Reviews.csv",
header = TRUE,
nrows = 8000)1 Pré-processamento
1.0.1 Bag-of-words
corp = VCorpus(VectorSource(dados$Text))
dtm = DocumentTermMatrix(corp,
control = list(tolower
= TRUE,
stemming
= FALSE,
removeNumbers
= TRUE,
removePunctuation
= TRUE,
removeStripwhitespace = TRUE,
weighting
= weightTf,
bounds=list(global=c(50, Inf))))
dtm## <<DocumentTermMatrix (documents: 8000, terms: 1045)>>
## Non-/sparse entries: 279916/8080084
## Sparsity : 97%
## Maximal term length: 13
## Weighting : term frequency (tf)dtm.matrix <- dtm %>% as.matrix # perdemos esparsidade :(2 Data-splitting
split <- sample(c("Treinamento",
"Validacao",
"Teste"),prob=c(0.6,0.2,0.2),size = nrow(dados),
replace = TRUE)3 Regressão
riscos <- list()3.1 Mínimos Quadrados
vc_lasso <- cv.glmnet(dtm.matrix[split=="Treinamento",], dados$Score[split=="Treinamento"],
alpha = 1)
predito_MMQ <- predict(vc_lasso, s = 0,
newx = dtm.matrix[split=="Teste",])Risco estimado:
riscos$minimos.quadrados <- (predito_MMQ-dados$Score[split=="Teste"])^2 %>% mean
riscos$minimos.quadrados## [1] 1.6074473.2 Lasso
plot(vc_lasso)
predito_lasso <- predict(vc_lasso, s = vc_lasso$lambda.min,
newx = dtm.matrix[split=="Teste",])Risco estimado:
riscos$lasso <- (predito_lasso-dados$Score[split=="Teste"])^2 %>% mean
riscos$lasso## [1] 1.286448Quais variáveis foram escolhidas?
table(coef(vc_lasso,s=vc_lasso$lambda.min)[,1]!=0)##
## FALSE TRUE
## 598 448coefs_estimates <- coef(vc_lasso,s=vc_lasso$lambda.min)
coefs <- coefs_estimates %>% as.matrix %>% as_tibble
names(coefs)="Estimativa"
coefs %>% mutate(variavel=rownames(coefs_estimates)) %>% arrange(desc(abs(Estimativa)))## # A tibble: 1,046 x 2
## Estimativa variavel
## <dbl> <chr>
## 1 4.07 (Intercept)
## 2 -1.18 worst
## 3 -1.10 terrible
## 4 -0.906 awful
## 5 -0.763 waste
## 6 -0.711 disappointed
## 7 -0.606 weak
## 8 -0.526 avoid
## 9 -0.516 return
## 10 -0.507 worse
## # ... with 1,036 more rowslibrary(ggpubr)
theme_set(theme_gray(base_size = 18))
coefs_estimates <- coef(vc_lasso,s=vc_lasso$lambda.min)
coefs_estimates = data.frame(Palavra = rownames(coefs_estimates), Coeficientes = coefs_estimates[,1])
coef_pos = coefs_estimates %>% arrange(desc(Coeficientes))
coef_neg = coefs_estimates %>% arrange(Coeficientes)
graf_pos = ggplot(data = coef_pos[2:30,],
aes(x = reorder(Palavra, Coeficientes), y = Coeficientes)) +
geom_bar(stat="identity", col = "white", fill = "dodgerblue1") + coord_flip() +
xlab("") + ylab("Estimativa") + ggtitle("Coeficientes Positivos")
graf_neg = ggplot(data = coef_neg[2:30,], aes(x = reorder(Palavra, -Coeficientes), y = Coeficientes)) +
geom_bar(stat="identity", col = "white", fill = "firebrick4") + coord_flip() +
xlab("Termo") + ylab("Estimativa") + ggtitle("Coeficientes Negativos")
ggarrange(graf_neg, graf_pos, ncol = 2, nrow = 1)
3.3 KNN
library(FNN)
k.grid <- round(seq(1,50,length.out = 10))
erro <- rep(NA,length(k.grid))
for(ii in seq_along(k.grid))
{
predito <- knn.reg(train=dtm.matrix[split=="Treinamento",],
y=dados$Score[split=="Treinamento"],
k=k.grid[ii])$pred
erro[ii] <- mean((predito-dados$Score[split=="Treinamento"])^2)
}
plot(k.grid,erro)
best.k <- k.grid[which.min(erro)]
predito <- knn.reg(train=dtm.matrix[split=="Treinamento",],
test = dtm.matrix[split=="Teste",],
y=dados$Score[split=="Treinamento"],
k=best.k)$predriscos$knn <- (predito-dados$Score[split=="Teste"])^2 %>% mean
riscos$knn## [1] 1.5647833.4 Árvores
library(rpart)
library(rpart.plot)
dados.tudo <- dtm.matrix %>% as.data.frame %>% mutate(Score=dados$Score)
# ajuste da arvore
arvore = rpart(Score ~.,
data = dados.tudo[split=="Treinamento",])
# poda
melhorCp = arvore$cptable[which.min(arvore$cptable[,"xerror"]),
"CP"]
poda = prune(arvore, cp = melhorCp)
# arvore
rpart.plot(poda, type = 4)
# predito árvore
predito_arvore = predict(poda,
dados.tudo[split=="Teste", ])riscos$arvores <- (predito_arvore-dados$Score[split=="Teste"])^2 %>% mean
riscos$arvores## [1] 1.5807993.5 Florestas
library(randomForest)
names(dados.tudo) <- make.names(names(dados.tudo))
# floresta aleatória
floresta = randomForest(Score ~.,
data = dados.tudo[split=="Treinamento",],
ntree = 100)
# predito floresta
predito_floresta = predict(floresta,
dados.tudo[split=="Teste", ])varImpPlot(floresta)
riscos$florestas <- (predito_floresta-dados$Score[split=="Teste"])^2 %>% mean
riscos$florestas## [1] 1.3223293.6 Comparando os métodos
library(knitr)
kable(as.matrix(unlist(riscos)),row.names = TRUE,digits=3)| minimos.quadrados | 1.607 |
| lasso | 1.286 |
| knn | 1.565 |
| arvores | 1.581 |
| florestas | 1.322 |
4 Classificação
riscos <- list()Definindo a variável resposta
dados$Bom <- ifelse(dados$Score>=4,"Bom",
"Ruim")4.1 Logística
vc_glm.lasso <- cv.glmnet(dtm.matrix[split=="Treinamento",], dados$Bom[split=="Treinamento"], alpha=1,family="binomial")
predito_logistica <- predict(vc_glm.lasso, s = 0,
newx = dtm.matrix[split=="Teste",],type="response")
hist(predito_logistica)
Avaliando o risco
classe_log <- ifelse(predito_logistica<1/2,"Bom","Ruim")
riscos$logistica <-mean(classe_log!=dados$Bom[split=="Teste"])
riscos$logistica## [1] 0.2292909table(classe_log,dados$Bom[split=="Teste"])##
## classe_log Bom Ruim
## Bom 959 172
## Ruim 174 2044.2 Linear
vc_lm.lasso <- cv.glmnet(dtm.matrix[split=="Treinamento",], dados$Bom[split=="Treinamento"]=="Bom", alpha=1)
predito_linear <- predict(vc_lm.lasso, s = 0,
newx = dtm.matrix[split=="Teste",],type="response")
hist(predito_linear)
Avaliando o risco
classe_lin <- ifelse(predito_linear>1/2,"Bom","Ruim")
riscos$linear <-mean(classe_lin!=dados$Bom[split=="Teste"])
riscos$linear## [1] 0.1868787table(classe_lin,dados$Bom[split=="Teste"])##
## classe_lin Bom Ruim
## Bom 1038 187
## Ruim 95 1894.3 Logística e linear com penalização
predito_logistica.p <- predict(vc_glm.lasso, s = vc_glm.lasso$lambda.min,
newx = dtm.matrix[split=="Teste",],type="response")
classe_log.p <- ifelse(predito_logistica.p<1/2,"Bom","Ruim")
riscos$logisitica.penalizada <-mean(classe_log.p!=dados$Bom[split=="Teste"])
riscos$logisitica.penalizada## [1] 0.1802518table(classe_log.p,dados$Bom[split=="Teste"])##
## classe_log.p Bom Ruim
## Bom 1072 211
## Ruim 61 165predito_linear.p <- predict(vc_lm.lasso, s = vc_lm.lasso$lambda.min,
newx = dtm.matrix[split=="Teste",],type="response")
classe_lin.p <- ifelse(predito_linear.p>1/2,"Bom","Ruim")
riscos$linear.penalizado <-mean(classe_lin.p!=dados$Bom[split=="Teste"])
riscos$linear.penalizado## [1] 0.1935056table(classe_lin.p,dados$Bom[split=="Teste"])##
## classe_lin.p Bom Ruim
## Bom 1092 251
## Ruim 41 1254.4 Naive Bayes
library(e1071)
n.b <- naiveBayes(x=dtm.matrix[split=="Treinamento",],y=as.factor(dados$Bom[split=="Treinamento"]))
predito_n.b <- predict(n.b,dtm.matrix[split=="Teste",])riscos$naive.bayes <-mean(predito_n.b!=dados$Bom[split=="Teste"])
riscos$naive.bayes## [1] 0.2485089table(predito_n.b,dados$Bom[split=="Teste"])##
## predito_n.b Bom Ruim
## Bom 919 161
## Ruim 214 2154.5 Árvores de Classificação
dados.tudo <- as.data.frame(dtm.matrix %>% as.data.frame %>% mutate(Bom=as.factor(dados$Bom)))
# ajuste da arvore
arvore = rpart(Bom ~., method = "class",
data = dados.tudo[split=="Treinamento",])
predito_arvore = predict(arvore, dados.tudo[split=="Teste",],type="class")
riscos$arvore <-mean(predito_arvore!=dados$Bom[split=="Teste"])
riscos$arvore## [1] 0.2491716table(predito_arvore,dados$Bom[split=="Teste"])##
## predito_arvore Bom Ruim
## Bom 1133 376
## Ruim 0 0Problema: classe muito desbalanceada!
mean(dados$Bom[split=="Teste"]=="Ruim")## [1] 0.2491716Alternativa: usar floresta de regressão para estimar \(P(\mbox{Bom}|x)\).
4.6 Florestas de Classificação
names(dados.tudo) <- make.names(names(dados.tudo))
# floresta aleatória
floresta = randomForest(Bom ~., method = "class",data = dados.tudo[split=="Treinamento",])
varImpPlot(floresta)
# predito floresta
predito_floresta = predict(floresta, dados.tudo[split=="Teste",])
table(predito_floresta)## predito_floresta
## Bom Ruim
## 1426 83riscos$floresta <-mean(predito_floresta!=dados$Bom[split=="Teste"])
riscos$floresta## [1] 0.2113983table(predito_floresta,dados$Bom[split=="Teste"])##
## predito_floresta Bom Ruim
## Bom 1120 306
## Ruim 13 704.7 KNN
k.grid <- round(seq(1,50,length.out = 10))
erro <- rep(NA,length(k.grid))
for(ii in seq_along(k.grid))
{
predito <- knn.cv(train=dtm.matrix[split=="Treinamento",], cl=dados$Bom[split=="Treinamento"],k=k.grid[ii])[1:sum(split=="Treinamento")]
erro[ii] <- mean((predito!=dados$Bom[split=="Treinamento"]))
}
plot(k.grid,erro)
best.k <- k.grid[which.min(erro)]
predito <- knn(train=dtm.matrix[split=="Treinamento",],cl=dados$Bom[split=="Treinamento"],
test = dtm.matrix[split=="Teste",],
k=best.k)riscos$knn <- mean(predito!=dados$Bom[split=="Teste"])
riscos$knn## [1] 0.2432074table(predito,dados$Bom[split=="Teste"])##
## predito Bom Ruim
## Bom 1115 349
## Ruim 18 274.8 XGBoost
set.seed(1)
dados <- fread("../../../../mineracaoDeDados/2016_01/notasDeAula/data/Reviews.csv",
header = TRUE,
nrows = 100000)
corp = VCorpus(VectorSource(dados$Text))
dtm = DocumentTermMatrix(corp,
control = list(tolower
= TRUE,
stemming
= FALSE,
removeNumbers
= TRUE,
removePunctuation
= TRUE,
removeStripwhitespace = TRUE,
weighting
= weightTf,
bounds=list(global=c(50, Inf))))
dtm## <<DocumentTermMatrix (documents: 100000, terms: 5640)>>
## Non-/sparse entries: 4392634/559607366
## Sparsity : 99%
## Maximal term length: 15
## Weighting : term frequency (tf)Proporção do teste menor!
split <- sample(c("Treinamento",
"Validacao",
"Teste"),prob=c(0.9,0.05,0.05),size = nrow(dados),
replace = TRUE)
table(split)## split
## Teste Treinamento Validacao
## 4997 89973 5030Converter para formato esparso aceito pelo xgboost.
library(Matrix)
nomes.variaveis=colnames(dtm)
dtm <- sparseMatrix(i=dtm$i, j=dtm$j, x=dtm$v, dims=c(dtm$nrow, dtm$ncol))dados$Bom <- as.numeric(dados$Score>=4)library(xgboost)
bstCV = xgb.cv(data = dtm[split=="Treinamento",],
label = dados$Bom[split=="Treinamento"],
nthread
= 7, nround
= 2000,
nfold
= 4, metrics = "error",
objective = "binary:logistic",
lambda
= .01, verbose = FALSE)
bst = xgboost(data = dtm[split=="Treinamento",],
label = dados$Bom[split=="Treinamento"], nthread = 7,
nround = which.min(bstCV$evaluation_log$test_error_mean),
objective = "binary:logistic", verbose = FALSE,
lambda = 0.01)
#saveRDS(bst,file = "bst.RDS")
#bst=loadRDS(file = "bst.RDS")importance_matrix=xgb.importance(nomes.variaveis, model = bst)[1:20,]
xgb.plot.importance(importance_matrix)
predicoes= predict(bst, dtm[split=="Teste",])
predito_bst = as.numeric(predicoes >1/2)
riscos$xgboost <- mean(predito_bst!=dados$Bom[split=="Teste"])
table(predito_bst,dados$Bom[split=="Teste"])##
## predito_bst 0 1
## 0 836 184
## 1 321 3656riscos$xgboost## [1] 0.1010606Aumentar \(n\) para a logística não é tão importante, mas para um método não paramétrico sim!
4.9 Comparando os métodos
library(knitr)
kable(as.matrix(unlist(riscos)),row.names = TRUE,digits=3)| logistica | 0.229 |
| linear | 0.187 |
| logisitica.penalizada | 0.180 |
| linear.penalizado | 0.194 |
| naive.bayes | 0.249 |
| arvore | 0.249 |
| floresta | 0.211 |
| knn | 0.243 |
| xgboost | 0.101 |