Muito se ouve falar que O Ensemble Learning ajuda a melhorar os resultados do aprendizado de máquina combinando vários modelos, e que essa abordagem permite a produção de melhor desempenho preditivo em comparação com um único modelo. Neste post, apresento o algoritmo Stacking, combinando alguns modelos interessantes como: Elastic NET, LASSO, Ridge e SVR, Spam.
Os Métodos de Ensemble são nada mais que meta-algoritmos que combinam várias técnicas de aprendizado de máquina em um único modelo preditivo para diminuir variação (bagging), viés (boosting) ou melhorar previsões (stacking).
Como esse post é voltado para melhora de previsões, os próximos passos são focados no método de Ensemble Learning via Stacking.
O Stacking é uma técnica de Ensemble Learning que combina vários modelos de classificação ou regressão por meio de um meta-classificador ou meta-regressor. Os modelos de nível básico são treinados com base em um conjunto de treinamento completo; o metamodelo é treinado nas saídas do modelo de nível básico como recursos. O nível básico geralmente consiste em diferentes algoritmos de aprendizado e, portanto, os conjuntos de empilhamento geralmente são heterogêneos.
Em outras palavras, basicamente o método Stacking consiste em tomar a amostra inteira, tomar rótulos via modelos de Machine Learning, sem combinação (ou seja, fazer previsões com esses modelos) e posteriormente, utilizando esses rótulos como entradas de um “novo conjunto de dados”, prever a variável resposta dos dados originais.
Talvez, o esquema visual mostrado abaixo possa te ajudar. Perceba que cada modelo simples (em vermelho) é aplicado individualmente aos dados orginais, a fim de prever a variável resposta (criar rótulos). Perceba também que esse rótulo estimado por cada modelo utilizado atuará como uma “variável explicativa” (entradas) para meu novo banco de dados. Por fim, o modelo combinado é apenas aplicar a um modelo qualquer (pode inclusive ser um dos mesmos utilizados individualmente) os dados formados pelos rótulos a fim de predizer o rótulo dos dados originais.
library(visNetwork)
modelos = data.frame(id = 1:7, group = c(rep("Modelo Simples",6),"Modelo Combinado"))
diagrama = data.frame(from = c(1:6), to = 7)
visNetwork(modelos, diagrama, width = "100%") %>%
visGroups(groupname = "Modelo Combinado", shape = "icon", icon = list(code = "f0c0", size = 75)) %>%
visGroups(groupname = "Modelo Simples", shape = "icon", icon = list(code = "f007", color = "red")) %>%
visLegend() %>%
addFontAwesome()Aqui, utilizarei um conjunto de dados simulados para comparar a performance da aplicação dos modelos de modo tradicional, contra o método Stacking. Os modelos escolhidos para ilustração dessa abordagem foram: SVR com Kernel Gaussiano, SpAM, Ridge, LASSO, Elastic Net e Floresta Aleatória. A métrica de desempenho utilizada foi a Estatística de Theil (1958). A Estatística proposta por Theil (1958) é uma derivação do estudo sobre o EQM, que pela sua fácil interpretabilidade, torna a comunicação entre corpo técnico e usuários com pouco expertise na área de ciência de dados, mais facilitada. Ahlburg, D. (1984), Pindyck, R. e Rubinfeld, D. (1997) e Makridakis, S.,Wheelwright, S. e Hyndman, R. (1998), são alguns trabalhos interessantes que tratam de aplicações utilizando a decomposição dessa métrica em três parcelas conhecidas como Proporções da Estatística de Theil, representadas da seguinte forma:
\[ U_m = \frac{(\bar{Y}- \bar{\hat{Y}})^2}{\frac{1}{n}\sum_{i=1}^n \left[ Y_i-\hat{Y}_i \right]^2} \ \ ; \ \ U_r = \frac{(\sigma_Y - \sigma_{\hat{Y}})^2}{\frac{1}{n}\sum_{i=1}^n \left[ Y_i-\hat{Y}_i \right]^2} \ \ ; \ \ U_c = \frac{2(1-r)\sigma_Y\sigma_{\hat{Y}}}{ \frac{1}{n}\sum_{i=1}^n \left[ Y_i-\hat{Y}_i \right]^2}\]
Perceba que por construção, \(U_m + U_r + U_c =1\), onde quanto mais \(U_c \approx 1\), melhor a performance do modelo, e consequentemente, menor o viés \(U_m\) e a varibilidade \(U_r\) dos resultados preditos. Ou seja, além da interpretação direta por meio de \(U_c\), a utilização das “Proporções da Estatística de Theil”, permitem a análise sobre o trade-off entre viés e variância.
A idéia foi construir um banco de dados pequeno, porém bastante heterogêneo, permitindo uma melhor análise sobre a diferença em performance entre os dois métodos.
n = 100
x1 = rnorm(n, 3,1)
x2 = rgamma(n,10,4)
x3 = runif(n,3,6)
x4 = rnorm(n,2,1)^2
x5 = log(runif(n,1,2))
y = x1 + x2 +x3 + x4 +x5
dados = data.frame(x1,x2,x3,x4,x5,y)Aqui, apresento funções auxiliares para aplicação dos modelos. A função medida trás a formulação das Proporções da Estatística de Theil.
varResp = 'y'
Medidas = function(residuos,pred,teste, formula){
ind = which(colnames(teste) == as.character(formula)[2])
y = teste[, ind]
estimados.medios= mean(pred)
observados.medios= mean(y)
estimados.desvios = sd(pred)
observados.desvios = sd(y)
r = cor(pred,y)
EQM = mean((residuos)^2)
numerador.Var = (estimados.medios - observados.medios)^2
numerador.Vies = (estimados.desvios - observados.desvios)^2
numerador.cov = 2*(1-r)*estimados.desvios*observados.desvios
Um = numerador.Vies/ EQM
Ur = numerador.Var/ EQM
Uc = numerador.cov/EQM
return(list(Um = Um,Ur = Ur,Uc = Uc))
}
formula = as.formula(paste(varResp, '~ .'))
AuxSaida = function(pred, teste, formula, alg){
ind = which(colnames(teste) == as.character(formula)[2])
residuos = pred-teste[, ind]
rs = Medidas(residuos,pred,teste,formula)
rs[['Classificador']] = alg
return(do.call(data.frame, rs))
}Aqui, apresento funções referentes ao processo de data splitting. Neste trabalho, fiz a validação baseando-se na validação cruzada k-fold repetida. Para um k=10, e 10 repetições.
ValidCruzada = function(cls, dados, formula, kfolds = 10, ...){
library(dismo)
id = kfold(1:nrow(dados), kfolds)
vcs = data.frame()
for(i in 1:kfolds){
treino = dados[id != i, ]
teste = dados[id == i, ]
kcls = cls(treino, teste, formula)
vcs = rbind(vcs, kcls)
}
vcs = vcs[, c(ncol(vcs), 1:(ncol(vcs) - 1))]
avg = aggregate(. ~ Classificador, data = vcs, FUN = mean)
std = aggregate(. ~ Classificador, data = vcs, FUN = sd)
return(list(Media = avg, Desvio = avg, Modelos = vcs))
}
ValidCruzadarep = function(cls, dados, formula, kfolds = 10, reps = 10){
x = data.frame()
for(i in 1:reps) x = rbind(x, ValidCruzada(cls, dados, formula, kfolds)$Media)
avg = aggregate(. ~ Classificador, data = x, FUN = mean)
std = aggregate(. ~ Classificador, data = x, FUN = sd)
return(list(Media = avg, Desvio = std, Modelos = x))
}#SVR
library(kernlab)
KERNEL2= "rbfdot"
SVR2 = function(treino, teste, formula){
cls = ksvm(formula, data = treino, kernel = KERNEL2)
pred = predict(cls, newdata = data.frame(teste))
return(pred)
}
medSVR2 = function(treino, teste, formula){
alg = 'SVR - rbfdot'
pred = SVR2(treino, teste, formula)
return(AuxSaida(pred, teste, formula, alg))
}
#SpAM
SPAM = function(treino, teste, formula){
library(SAM)
ind.treino = which(colnames(treino) == as.character(formula)[2])
ind.teste = which(colnames(teste) == as.character(formula)[2])
dados.x.treino= (treino[,-ind.treino])
dados.y.treino= (treino[,ind.treino])
dados.x.teste = (teste[,-ind.teste])
for(i in 1:ncol(dados.x.treino)) {
dados.x.treino[,i] <- as.numeric(dados.x.treino[,i])
}
for(i in 1:ncol(dados.x.teste)) {
dados.x.teste[,i] <- as.numeric(dados.x.teste[,i])
}
cls = samQL(dados.x.treino,dados.y.treino,p=1)
pred = predict(object=cls, dados.x.teste)
pred= matrix(pred$values,nrow= length(teste[,ind.teste]),
ncol=(length(pred$values)/length(teste[,ind.teste])))
pred = apply(pred,1, mean)
return(pred)
}
medSPAM = function(treino, teste, formula){
alg = 'SpAM'
pred = SPAM(treino, teste, formula)
return(AuxSaida(pred, teste, formula, alg))
}
#Ridge
library(glmnet)
RIDGE= function(treino, teste, formula){
ind.treino = which(colnames(treino) == as.character(formula)[2])
ind.teste = which(colnames(teste) == as.character(formula)[2])
dados.x.treino= data.matrix(treino[,-ind.treino])
dados.y.treino= data.matrix(treino[,ind.treino])
dados.x.teste = data.matrix(teste[,-ind.teste])
cls = cv.glmnet(dados.x.treino,dados.y.treino,alpha = 0,type.measure = "mse")
lambda= as.numeric(cls$lambda.1se)
pred = predict(object=cls, s =lambda, newx =dados.x.teste)
return(pred)
}
medRIDGE = function(treino, teste, formula){
alg = 'RIDGE'
pred = RIDGE(treino, teste, formula)
return(AuxSaida(pred, teste, formula, alg))
}
#LASSO
library(glmnet)
LASSO= function(treino, teste, formula){
ind.treino = which(colnames(treino) == as.character(formula)[2])
ind.teste = which(colnames(teste) == as.character(formula)[2])
dados.x.treino= data.matrix(treino[,-ind.treino])
dados.y.treino= data.matrix(treino[,ind.treino])
dados.x.teste = data.matrix(teste[,-ind.teste])
cls = cv.glmnet(dados.x.treino,dados.y.treino,alpha = 1,type.measure = "mse")
lambda= as.numeric(cls$lambda.1se)
pred = predict(object=cls, s =lambda, newx =dados.x.teste)
return(pred)
}
medLASSO = function(treino, teste, formula){
alg = 'LASSO'
pred = LASSO(treino, teste, formula)
return(AuxSaida(pred, teste, formula, alg))
}
#ELASTIC NET
library(glmnet)
E.NET= function(treino, teste, formula){
ind.treino = which(colnames(treino) == as.character(formula)[2])
ind.teste = which(colnames(teste) == as.character(formula)[2])
dados.x.treino= data.matrix(treino[,-ind.treino])
dados.y.treino= data.matrix(treino[,ind.treino])
dados.x.teste = data.matrix(teste[,-ind.teste])
cls = cv.glmnet(dados.x.treino,dados.y.treino,alpha = 0.5,type.measure = "mse")
lambda= as.numeric(cls$lambda.1se)
pred = predict(object=cls, s =lambda, newx =dados.x.teste)
return(pred)
}
medE.NET = function(treino, teste, formula){
alg = 'ELASTIC NET'
pred = E.NET(treino, teste, formula)
return(AuxSaida(pred, teste, formula, alg))
}
#FLORESTAS ALEATæ¼ã¸³RIAS
Floresta = function(treino, teste, formula){
library(randomForest)
cls = randomForest(formula, data = data.frame(treino),importance= TRUE,
proximity=TRUE)
pred = predict(cls, newdata = data.frame(teste))
return(pred)
}
medFloresta = function(treino, teste, formula){
alg = 'Floresta Aleatæ¼ã¸³ria'
pred = Floresta(data.frame(treino), data.frame(teste), formula)
return(AuxSaida(pred, teste, formula, alg))
}tabela.de.desempenho = matrix(NA, ncol = 4, nrow = 6)
colnames(tabela.de.desempenho) = c("Modelos","Um","Ur", "Uc")
tabela.de.desempenho = data.frame(tabela.de.desempenho)
medida = medSVR2
nome.medida = "SVR - Kernel Gaussiano"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho[1,] = c(nome.medida,round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medE.NET
nome.medida = "Elastic Net"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho[2,] = c(nome.medida, round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medLASSO
nome.medida = "Lasso"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho[3,] = c(nome.medida,round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medRIDGE
nome.medida = "Ridge"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho[4,] = c(nome.medida,round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medSPAM
nome.medida = "SpAM"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho[5,] = c(nome.medida,round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medFloresta
nome.medida = "Floresta Aleatæ¼ã¸³ria"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho[6,] = c(nome.medida, round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
tabela.de.desempenho= tabela.de.desempenho[order(tabela.de.desempenho[,4]),]
library(formattable)
tabela.de.desempenho = data.frame(tabela.de.desempenho)
formattable(tabela.de.desempenho, list(Uc =color_tile("white", "red")))| Modelos | Um | Ur | Uc | |
|---|---|---|---|---|
| 4 | Ridge | 0.945 | 0.119 | 0.033 |
| 5 | SpAM | 0.754 | 0.144 | 0.197 |
| 6 |
Floresta Aleat |
0.695 | 0.107 | 0.296 |
| 2 | Elastic Net | 0.671 | 0.122 | 0.305 |
| 3 | Lasso | 0.495 | 0.141 | 0.46 |
| 1 | SVR - Kernel Gaussiano | 0.333 | 0.1 | 0.667 |
Essa seção está dividia em duas etapas: na primeira construímos o banco de dado, conforme descrito na Seção 1; por fim, na segunda etapa, fazemos tunagem sobre os modelos estudados a fim de encontrar o melhor regressor.
library(knitr)
library(kableExtra)
treino = dados
teste = dados
###########################################################
y_floresta = Floresta(treino, teste, formula)
y_E.NET = c(E.NET(treino, teste, formula))
y_LASSO = c(LASSO(treino, teste, formula))
y_RIDGE = c(RIDGE(treino, teste, formula))
y_SPAM = SPAM(treino, teste, formula)
y_SVR2 = SVR2(treino, teste, formula)
dados.preditos = data.frame('Floresta' = y_floresta, 'E.Net' = y_E.NET,
'Lasso' = y_LASSO, 'Ridge' = y_RIDGE, 'Spam' = y_SPAM,
"SVR" = y_SVR2,"Y" = dados$y)
kable(dados.preditos[1:10,], format = "markdown", digits = 2)| Floresta | E.Net | Lasso | Ridge | Spam | SVR | Y |
|---|---|---|---|---|---|---|
| 14.36 | 14.51 | 14.50 | 14.51 | 14.46 | 14.17 | 14.55 |
| 15.68 | 16.12 | 16.08 | 16.05 | 14.91 | 16.67 | 16.29 |
| 20.63 | 23.35 | 23.40 | 22.73 | 21.92 | 23.06 | 23.51 |
| 12.01 | 12.05 | 12.04 | 12.17 | 12.43 | 12.09 | 12.05 |
| 9.95 | 8.10 | 8.12 | 8.52 | 9.99 | 8.90 | 7.86 |
| 19.60 | 20.93 | 20.99 | 20.46 | 19.92 | 21.21 | 20.99 |
| 13.73 | 14.03 | 14.03 | 14.03 | 14.23 | 13.63 | 14.01 |
| 12.49 | 12.56 | 12.53 | 12.70 | 12.73 | 12.26 | 12.58 |
| 18.77 | 19.94 | 19.91 | 19.66 | 18.27 | 19.51 | 20.16 |
| 13.35 | 13.51 | 13.53 | 13.54 | 13.77 | 13.30 | 13.45 |
Quando não sabemos qual modelo atua melhor como meta-regressor é interessagem testar uma série de modelos escolhendo o de melhor performance. Neste caso, fiz uso dos mesmo modelos aplicados anteriormente. Isto facilita a comparação entre os métodos tradicional e stacking, uma vez que estaremos comparando os mesmos modelos, mas aplicados à diferentes algoritmos.
dados = dados.preditos
varResp = 'Y'
formula = as.formula(paste(varResp, '~ .'))
tabela.de.desempenho.stacking = matrix(NA, ncol = 4, nrow = 6)
colnames(tabela.de.desempenho.stacking) = c("Modelos","Um","Ur", "Uc")
tabela.de.desempenho.stacking = data.frame(tabela.de.desempenho.stacking)
medida = medSVR2
nome.medida = "SVR - Kernel Gaussiano"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho.stacking[1,] = c(nome.medida,round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medE.NET
nome.medida = "Elastic Net"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho.stacking[2,] = c(nome.medida, round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medLASSO
nome.medida = "Lasso"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho.stacking[3,] = c(nome.medida,round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medRIDGE
nome.medida = "Ridge"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho.stacking[4,] = c(nome.medida,round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medSPAM
nome.medida = "SpAM"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho.stacking[5,] = c(nome.medida,round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
medida = medFloresta
nome.medida = "Floresta Aleatæ¼ã¸³ria"
saida = ValidCruzadarep(medida, dados, formula, kfolds = 10, reps = 10)
tabela.de.desempenho.stacking[6,] = c(nome.medida, round(saida$Media[[2]],3),
round(saida$Media[[3]],3),round(saida$Media[[4]],3))
tabela.de.desempenho.stacking= tabela.de.desempenho.stacking[order(tabela.de.desempenho.stacking[,4]),]
tabela.de.desempenho.stacking = data.frame(tabela.de.desempenho.stacking)
formattable(tabela.de.desempenho.stacking, list(Uc =color_tile("white", "red")))| Modelos | Um | Ur | Uc | |
|---|---|---|---|---|
| 5 | SpAM | 0.979 | 0.109 | 0.011 |
| 3 | Lasso | 0.733 | 0.111 | 0.255 |
| 1 | SVR - Kernel Gaussiano | 0.436 | 0.169 | 0.488 |
| 4 | Ridge | 0.283 | 0.114 | 0.702 |
| 2 | Elastic Net | 0.286 | 0.104 | 0.71 |
| 6 |
Floresta Aleat |
0.24 | 0.121 | 0.737 |
Perceba que podemos destacar um ganho preditivo, gerado principalmente pela redução da parcela da Estatistica de Theil destinada ao viés, entre a melhor aplicação tradicional (Modelo SVR com Kernel Gaussiano) e o melhor método de aprendizagem via Stacking (pelo modelo Floresta Aleatória).
Vale ressaltar que nem sempre o ganho é significativo o sufiente para que o método Stacking seja utilizado. Isto porque, além de elevar o custo computacional, pode resultar em problemas como overffintg. Afinal, aplicar Stacking é elevar a complexidade do modelo final, e quando isso é feito sem uma análise criteriosa pode tornar seu modelo pouco generalista.