Como Fazer Stacking de Modelos de Machine Learning

Quem acompanha competições sabe que uma das coisas mais importantes é saber juntar vários modelos para criar uma solução poderosa. Várias pessoas já me perguntaram, por e-mail ou nas apresentações que fiz, sobre ensembles. Este é um assunto importante não apenas para competições, mas também para casos reais onde se quer extrair o máximo possível de performance dos modelos.

Ensembles são conjuntos de modelos que oferecem uma performance melhor do que cada modelo que o compõe.

Então neste artigo quero exemplificar a melhor maneira que conheço de criar ensembles: stacking. Este é um método que usei em todas as competições que tive bons resultados.

Antes de começar, uma dica: é importante pensar em Machine Learning como “processos”. Aqui não vamos apenas testar modelos em um banco de dados, mas vamos testar processos, pipelines de métodos aplicadas aos dados para saber quais são os resultados.

Este material é parte da apresentação que fiz no primeiro PyData Meetup São Paulo. O Jupyter Notebook original está diponível neste link.

import pandas as pd
import numpy as np

Índice

• Carregando os dados
• Conjuntos de Variáveis (Features)
  • Funções auxiliares
  • Feature set 1: variáveis “numéricas”
  • Feature set 2: Ordinal Encoding Categóricas
  • Bônus: OneHot Encoding Categóricas
• Transformações do Target
  • Log
  • Raiz Quadrada
• Gerando modelos com modelos/algoritmos diferentes
  • Ajustando hiperparâmetros
• Stacking

Carregando os dados

Os dados utilizados são transações comerciais de imóveis na cidade de Ames, Iowa. Nosso objetivo é prever o preço de venda de uma casa alimentando o modelo com as características. Estes dados também são tema de uma competição no Kaggle. Os dados de treino e teste podem ser encontrados em: https://www.kaggle.com/c/house-prices-advanced-regression-techniques/

Algo muito importante é explorar os dados em busca de ideias de features e métodos para validar os modelos. Como o foco deste material é stacking, vou pular esta parte. Ainda assim, acho importante esclarecer que temos dois tipos básicos de variáveis nestes dados: categóricas e numéricas.

As variáveis categóricas possuem níveis que não podem ser ordenados. Em alguns casos é possível pensar em maneiras de ordená-los, como uma variável que descreve se uma rua é asfaltada ou não. Neste caso pode-se pensar que a rua asfaltada vai valorizar o imóvel.

As variáveis numéricas podem receber qualquer valor contínuo. Um exemplo seria o tamanho do terreno.

Nesta célula carregamos os dados e criamos um DataFrame com as features (X) e uma Série com o preço da venda (SalePrice, y) que é nosso alvo.

No total temos 79 colunas, mas limitei a 5 aqui por causa da formatação do site.

train = pd.read_csv('train.csv', index_col='Id')
X, y = train.drop('SalePrice', axis=1), train.SalePrice.copy()
train.head().iloc[:, :5]

Conjuntos de Variáveis (Features)

Uma das maneiras de obter modelos diferentes é variar a representação dos dados utilizada para treiná-los. Por isso nosso primeiro passo será construir estas representações.

Funções auxiliares

A métrica sugerida pelo Kaggle para avaliar os modelos é o RMSLE, este erro leva em conta a diferença entre o logaritmo das previsões e do alvo. É possível pensar neste erro como uma aproximação do erro percentual do modelo, mas com propriedades mais interessantes do ponto de vista matemático.

Criei três funções para cálculo do erro. Nós vamos fazer transformações da variável y, então para facilitar nosso trabalho, decidi criar funções personalizadas para cada transformação, assim evita alguma confusão na hora de transformá-las para computar o erro numa função só.

Para usarmos algumas funções do Scikit-learn que vão ajudar a manter nosso código mais limpo e mais eficiente, precisamos criar as funções de erro da maneira requerida pelo módulo. Neste caso, a função pecisa receber um modelo treinado, as features e o alvo. Ela computará as previsões e deverá retornar um número, que é o valor da métrica de erro.

A última linha da célula cria um objeto da validação cruzada do scikit-learn. Ele vai cuidar da divisão dos dados para nós. Este assunto merece uma série de artigos, mas basicamente é uma maneira de dividir os dados em N partes, e usar repetidamente N-1 partes como dados de treino, e a parte restante como dados de teste. Por exemplo, se temos as partes [1, 2, 3, 4, 5], numa das iterações podemos treinar usando as partes [1, 2, 3, 4] e validar na parte 5.

Apesar desta ser uma série temporal, o Kaggle decidiu tratar como dados aleatórios independentes, então por isso faremos esta validação cruzada simples.

from sklearn.metrics import mean_squared_error
from sklearn.model_selection import cross_val_score, cross_val_predict, KFold
from sklearn.ensemble import RandomForestRegressor

def rmsle(estimator, X, y):
    p = estimator.predict(X)
    return np.sqrt(mean_squared_error(np.log1p(y), np.log1p(p)))

def rmsle_log_y(estimator, X, y):
    p = estimator.predict(X)
    return np.sqrt(mean_squared_error(y, p))

def rmsle_sqrt_y(estimator, X, y):
    p = estimator.predict(X)
    y = np.power(y, 2)
    p = np.power(p, 2)
    return np.sqrt(mean_squared_error(np.log1p(y), np.log1p(p)))

kf = KFold(n_splits=5, shuffle=True, random_state=1)

Feature set 1: variáveis “numéricas”

O primeiro conjunto de features que teremos será uma seleção simples das variáveis originalmente numéricas dos dados. Para isto basta selecionar todas as colunas que possuem variáveis com números inteiros ou de ponto flutuante. Além disso, para substituir os valores nulos (o scikit-learn exige a substituição), decidi colocar o valor -1. Como vamos utilizar apenas modelos baseados em árvores com estes dados originais, não há necessidade de se preocupar muito com isso para nossos propósitos.

Após armazenar estes dados na variável X1, vemos que são 36 colunas. Já nesta parte quero treinar um modelo de Random Forest dentro da validação cruzada, para sabermos como ele se sairia sozinho nos dados originais. Para isso utilizei a função cross_val_score. Ela é uma função que facilita bastante na hora de fazer a validação cruzada com scikit-learn. Basta colocar um modelo e os dados. Neste caso optei por especificar um esquema de validação e uma métrica de erro personalizada.

Esta função retorna uma lista com os erros de cada iteração da validação cruzada, então fiz a média para sabermos o erro médio das partes.

X1 = X.select_dtypes(include=[np.number]).fillna(-1)
print('Dims', X1.shape)
model = RandomForestRegressor(n_estimators=1000, random_state=0)
error = cross_val_score(model, X1, y, cv=kf, scoring=rmsle).mean()
print('RMSLE:', error)

Dims (1460, 36)
RMSLE: 0.14582352618

Feature set 2: Ordinal Encoding Categóricas

Agora vamos criar um outro conjunto de features, desta vez adicionando as variáveis categóricas. Existem várias maneiras de codificar este tipo de variável para os modelos, uma delas é usando um formato ordinal. Isso simplesmente significa substituir cada valor original por números sequenciais. Em alguns modelos isso pode ser problemático, pois eles tentarão capturar alguma relação de ordem em valores que podem não ter. No nosso caso, com modelos baseados em árvores de decisão, este problema é quase inexistente.

Após codificar desta maneira, rodamos novamente a validação cruzada, agora nestes novos dados.

from sklearn.preprocessing import LabelEncoder

X2 = X.copy()
for col in X2.columns:
    if X2[col].dtype == object:
        enc = LabelEncoder()
        X2[col] = enc.fit_transform(X[col].fillna('Missing'))

print('Dims', X2.shape)
X2.fillna(-1, inplace=True)
model = RandomForestRegressor(n_estimators=1000, random_state=0)
error = cross_val_score(model, X2, y, cv=kf, scoring=rmsle).mean()
print('RMSLE:', error)

Dims (1460, 79)
RMSLE: 0.143837364859

Bônus: OneHot Encoding Categóricas

A maneira mais popular de codificar variáveis categóricas é o One Hot Encoding. Basicamente consiste em transformar cada valor da variável em uma coluna cujo novo valor será 1 caso a variável tenha aquele valor num determinado exemplo, ou 0 em caso negativo. Existem indicativos que árvores de decisão não processem tão bem este tipo de representação, mas em alguns casos práticos já vi isto funcionar melhor que o ordinal, então veja esta como mais uma ferramenta.

Este método cria mais de 200 novas colunas, o que deixa o processo de treino mais devagar, então decidi deixar a linha da validação cruzada comentada. Caso queira ver o resultado, basta rodá-la sem o jogo da velha.

#from sklearn.preprocessing import OneHotEncoder
X3 = X.copy()
cats = []
for col in X3.columns:
    if X3[col].dtype == object:
        X3 = X3.join(pd.get_dummies(X3[col], prefix=col), how='left')
        X3.drop(col, axis=1, inplace=True)
    

print('Dims', X3.shape)
X3.fillna(-1, inplace=True)
model = RandomForestRegressor(n_estimators=1000, random_state=0)
#cross_val_score(model, X3, y, cv=kf, scoring=rmsle).mean()

Dims (1460, 288)

Transformações do Target

Uma maneira interessante de criar diversidade, e às vezes até obter uma melhor performance, num caso de regressão, é transformar a variável que estamos tentando prever. Neste caso testaremos duas transformações: logaritmo e raiz quadrada.

Log

É possível ver que tentar prever o logaritmo do preço nos dá um resultado melhor. Isto acontece não só pelo fato do modelo capturar padrões diferentes, mas também porque usamos uma métrica baseada na diferença de logaritmos.

model = RandomForestRegressor(n_estimators=1000, random_state=0)
error = cross_val_score(model, X1, np.log1p(y), cv=kf, scoring=rmsle_log_y).mean()
print('RF, X1, log-target RMSLE:', error)

model = RandomForestRegressor(n_estimators=1000, random_state=0)
error = cross_val_score(model, X2, np.log1p(y), cv=kf, scoring=rmsle_log_y).mean()
print('RF, X2, log-target RMSLE:', error)

RF, X1, log-target RMSLE: 0.14518580749
RF, X2, log-target RMSLE: 0.14207134495

Raiz Quadrada

Esta transformação também nos dá um resultado melhor do que usar a variável em seu estado original. Uma das sugestões da razão pela qual vemos este efeito é que estas transformações fazem com que a variável y tenha uma distribuição mais próxima da normal, o que facilita o trabalho do modelo.

model = RandomForestRegressor(n_estimators=1000, random_state=0)
error = cross_val_score(model, X1, np.sqrt(y), cv=kf, scoring=rmsle_sqrt_y).mean()
print('RF, X1, sqrt-target RMSLE:', error)

model = RandomForestRegressor(n_estimators=1000, random_state=0)
error = cross_val_score(model, X2, np.sqrt(y), cv=kf, scoring=rmsle_sqrt_y).mean()
print('RF, X2, sqrt-target RMSLE:', error)

RF, X1, sqrt-target RMSLE: 0.145652934484
RF, X2, sqrt-target RMSLE: 0.143004600132

Gerando modelos com modelos/algoritmos diferentes

Outra maneira de gerar diversidade para o ensemble é gerar modelos diferentes. Neste caso vou usar meu modelo preferido, o GBM. Este também é baseado em árvores de decisão, mas basicamente treina cada árvore sequencialmente focando nos erros cometidos pelas anteriores.

Nas células abaixo é possível ver a performance deste modelo nos feature sets e transformações que usamos com a Random Forest. Vemos que ele traz uma melhora significativa, capturando melhor os padrões da relação entre as variáveis e o preço de venda dos imóveis.

from sklearn.ensemble import GradientBoostingRegressor
                
model = GradientBoostingRegressor(random_state=0)
error = cross_val_score(model, X1, np.log1p(y), cv=kf, scoring=rmsle_log_y).mean()
print('GBM, X1, log-target RMSLE:', error)

model = GradientBoostingRegressor(random_state=0)
error = cross_val_score(model, X2, np.log1p(y), cv=kf, scoring=rmsle_log_y).mean()
print('GBM, X2, log-target RMSLE:', error)

GBM, X1, log-target RMSLE: 0.133492454914
GBM, X2, log-target RMSLE: 0.129806890482

from sklearn.ensemble import GradientBoostingRegressor
                
model = GradientBoostingRegressor(random_state=0)
error = cross_val_score(model, X1, np.sqrt(y), cv=kf, scoring=rmsle_sqrt_y).mean()
print('GBM, X1, sqrt-target RMSLE:', error)

model = GradientBoostingRegressor(random_state=0)
error = cross_val_score(model, X2, np.sqrt(y), cv=kf, scoring=rmsle_sqrt_y).mean()
print('GBM, X2, sqrt-target RMSLE:', error)

GBM, X1, sqrt-target RMSLE: 0.134258972813
GBM, X2, sqrt-target RMSLE: 0.130919235682

Ajustando hiperparâmetros

Para simplificar o exemplo e focar na parte do ensemble, não vou fazer o ajuste dos hiperparâmetros do modelo. Hiperparâmetros são os atributos do modelo (como a profundidade das árvores de decisão) que precisam ser ajustados usando dados separados de validação ou um ciclo de validação cruzada.

É bom saber que nem sempre os melhores modelos formam os melhores ensembles. É importante ter modelos poderosos ao fazer stacking, mas também devemos lembrar que é importante ter diversidade. Às vezes alguns modelos que possuem um erro mais alto podem capturar padrões diferentes dos melhores modelos e por isso contribuem com o ensemble.

Caso você decida fazer o ajuste dos hiperparâmetros, é importante colocá-lo dentro do ciclo de validação cruzada que veremos no próximo passo.

Stacking

Tudo o que fizemos acima é para que possamos criar o nosso ensemble. Esta é a hora juntarmos os métodos usados para melhorarmos o poder preditivo dos nossos modelos.

O Stacking é uma maneira de fazer o ensemble na qual usamos modelos para fazer previsões, e depois usamos estas previsões como features em novos modelos, no que pode ser chamado de “segundo nível”. Você pode fazer este processo várias vezes, mas a cada nível o retorno em performance com relação à computação necessária é menor.

Nesta fase precisamos de dois ciclos de validação cruzada: externo e interno. No interno, nós treinaremos os modelos nos dados originais e faremos as previsões. No externo, treinaremos o modelo usando as previsões do primeiro passo como features.

A cada passo da validação cruzada interna, vamos salvar as previsões para a parte dos dados que for usada como validação. Desta maneira teremos previsões para todos os exemplos de nossos dados de treino. Além disso treinaremos um modelo nos dados originais de treino para podermos fazer previsões para os dados de teste.

No ciclo externo treinaremos um modelo nas previsões geradas pelo ciclo interno, e as features dos dados de validação serão as previsões dos modelos de primeiro nível nos dados de teste.

No nosso caso específico, criaremos previsões usando todas as combinações de modelos (RF e GBM), transformações do target (log e raiz quadrada) e feature sets (X1 e X2). Não usaremos o X3 porque levaria muito tempo para treinar, e para os fins de demonstração do método estes dois serão o bastante.

No fim teremos previsões de 8 modelos no primeiro nível. No segundo nível usei uma regressão linear regularizada (Ridge). Tendo estas previsões podemos computar o erro nos dados fora das nossas amostras de treino e validação internas, o que nos dará uma estimativa confiável do erro do ensemble.

from itertools import product
from sklearn.linear_model import Ridge

kf_out = KFold(n_splits=5, shuffle=True, random_state=1)
kf_in = KFold(n_splits=5, shuffle=True, random_state=2)

cv_mean = []
for fold, (tr, ts) in enumerate(kf_out.split(X, y)):
    X1_train, X1_test = X1.iloc[tr], X1.iloc[ts]
    X2_train, X2_test = X2.iloc[tr], X2.iloc[ts]
    y_train, y_test = y.iloc[tr], y.iloc[ts]
    
    modelos = [GradientBoostingRegressor(random_state=0), RandomForestRegressor(random_state=0)]
    targets = [np.log1p, np.sqrt]
    feature_sets = [(X1_train, X1_test), (X2_train, X2_test)]
    
    
    predictions_cv = []
    predictions_test = []
    for model, target, feature_set in product(modelos, targets, feature_sets):
        predictions_cv.append(cross_val_predict(model, feature_set[0], target(y_train), cv=kf_in).reshape(-1,1))
        model.fit(feature_set[0], target(y_train))
        ptest = model.predict(feature_set[1])
        predictions_test.append(ptest.reshape(-1,1))
    
    predictions_cv = np.concatenate(predictions_cv, axis=1)
    predictions_test = np.concatenate(predictions_test, axis=1)
    
    stacker = Ridge()
    stacker.fit(predictions_cv, np.log1p(y_train))
    error = rmsle_log_y(stacker, predictions_test, np.log1p(y_test))
    cv_mean.append(error)
    print('RMSLE Fold %d - RMSLE %.4f' % (fold, error))
    
print('RMSLE CV5 %.4f' % np.mean(cv_mean))
    

RMSLE Fold 0 - RMSLE 0.1248
RMSLE Fold 1 - RMSLE 0.1449
RMSLE Fold 2 - RMSLE 0.1257
RMSLE Fold 3 - RMSLE 0.1409
RMSLE Fold 4 - RMSLE 0.1087
RMSLE CV5 0.1290

Como podemos ver, nosso melhor modelo de primeiro nível é o GBM treinado com a transformação log nos dados X2, que atinge o erro de 0,1298. Nosso ensemble atinge o valor de 0,1290. Uma melhora de 0,62%.

O objetivo deste artigo era demonstrar o método, sem se preocupar muito com a performance no fim. Um ensemble feito visando melhora de performance pode apresentar um resultado mais significativo.

Em alguns casos, como em fundos de investimento ou aplicações de saúde, uma melhora pequena pode ter um resultado bastante significativo no mundo real, que justifique a criação de uma solução mais completa usando stacking.

Como sempre na aplicação de Machine Learning, nada é garantia de sucesso, mas este método é um dos mais consistentes em oferecer uma melhora.