Implementação de GLM com Grid Search no R usando o H2O.ai como backend

Para quem usa o R não existe nada mais irritante do que ter que lidar com o péssimo gerenciamento de memória da ferramenta, o que limita e muito o uso do R como uma ferramenta séria para a construção de modelos que possam ir para produção e possam permitir a construção de plataformas/sistemas inteligentes.

Vamos aqui em algumas linhas mostrar como usar o H2O.ai como backend de processamento (o que abstraí todos esses problemas de memória e processamento) para a criação de um modelo usando GLM.

O pulo do gato aqui é que o H2O faz todo o gerenciamento de memória, e independente da sua fonte de dados ele faz todo o pipeline do buffer de memória de forma que não há estouro de memória; ou mesmo uma lentidão generalizada no sistema.

Esse exemplo é baseado totalmente na documentação do H2O e tem o objetivo somente de mostrar como essa ferramenta funciona.

Nesse caso eu vou usar em um notebook, mas poderia ser utilizado por exemplo em uma máquina na Amazon usando o comando abaixo no momento da inicialização do cluster:

#Production Cluster (Not applicable)
#localH2O <- h2o.init(ip = '10.112.81.210', port =54321, nthreads=-1) # Máquina 1
#localH2O <- h2o.init(ip = '10.112.80.74', port =54321, nthreads=-1) # Máquina 2 - Sim, aqui usamos um cluster com dois nós para processamento! 😉

Primeiramente vamos remover qualquer instalação antiga do H2O.ai da máquina em questão:

# The following two commands remove any previously installed H2O packages for R.
if ("package:h2o" %in% search()) { detach("package:h2o", unload=TRUE) }
if ("h2o" %in% rownames(installed.packages())) { remove.packages("h2o") }

Em seguida vamos fazer o download e instalação de todos os pacotes dos quais o H2O tem alguma dependência direta ou indireta.

# Next, we download packages that H2O depends on.
if (! ("methods" %in% rownames(installed.packages()))) { install.packages("methods") }
if (! ("statmod" %in% rownames(installed.packages()))) { install.packages("statmod") }
if (! ("stats" %in% rownames(installed.packages()))) { install.packages("stats") }
if (! ("graphics" %in% rownames(installed.packages()))) { install.packages("graphics") }
if (! ("RCurl" %in% rownames(installed.packages()))) { install.packages("RCurl") }
if (! ("jsonlite" %in% rownames(installed.packages()))) { install.packages("jsonlite") }
if (! ("tools" %in% rownames(installed.packages()))) { install.packages("tools") }
if (! ("utils" %in% rownames(installed.packages()))) { install.packages("utils") }

Em seguida faremos a instalação da lib do H2O e o instanciamento da lib no R Studio.

# Now we download, install and initialize the H2O package for R.
install.packages("h2o", type="source", repos=(c("http://h2o-release.s3.amazonaws.com/h2o/rel-turing/8/R")))
# Load library
library(h2o)

Instalação feita e biblioteca carregada, vamos agora para algumas configurações.

No próprio R Studio você pode escolher o número de processadores no qual o cluster (nesse caso o seu notebook/desktop) vai utilizar. Lembrando que quanto maior for o número de cores utilizados, mais processamento o H2O vai consumir e menos recursos estarão disponíveis para outras tarefas. O padrão é a utilização de 2 cores, mas no meu caso eu vou usar todos os processadores.

# Start instance with all cores. 
# The -1 is the parameter to use with all cores. Use this carefully.
# The default parameter is 2 cores. 
h2o.init(nthreads = -1)

Agora vamos ver as informações do nosso cluster:

# Cluster Info
h2o.clusterInfo()

# R is connected to the H2O cluster: 
#   H2O cluster uptime:         3 seconds 267 milliseconds 
# H2O cluster version:        3.10.0.8 
# H2O cluster version age:    2 months and 26 days  
# H2O cluster name:           H2O_started_from_R_flavio.clesio_udy929 
# H2O cluster total nodes:    1 
# H2O cluster total memory:   1.78 GB 
# H2O cluster total cores:    4 
# H2O cluster allowed cores:  4 
# H2O cluster healthy:        TRUE 
# H2O Connection ip:          localhost 
# H2O Connection port:        54321 
# H2O Connection proxy:       NA 
# R Version:                  R version 3.3.2 (2016-10-31) 

Como podemos ver dos 4 processadores no total, estou usando todos eles (allowed cores) para o processamento.

Outro fato que podemos ver aqui é o que o H2O também está instanciado para usar a GUI na Web. Para isso, basta entrar no endereço no navegador com o endereço http://localhost:54321/flow/index.html.

Para este exemplo, vamos usar a base de dados Airlines que contém diversas informações reais de voos nos EUA e todas as causas de atraso de 1987 até 2008. A versão completa com 12Gb pode ser encontrada aqui.

Seguindo adiante, vamos agora fazer o carregamento dos dados direto de uma URL e importar em um objeto do R.

# GLM Demo Deep Dive
# Path of normalized archive. Can be a URL or a local path 
airlinesURL = "https://s3.amazonaws.com/h2o-airlines-unpacked/allyears2k.csv"
# We'll create the object .hex (extention of data files in H2O) 
# and using the importFile property, we'll set the path and the destination frame.
# As default behaviour H2O parse the file automatically.
airlines.hex = h2o.importFile(path = airlinesURL, destination_frame = "airlines.hex")

Neste caso o objeto airlines.hex é será o dataframe no qual o H2O irá aplicar os algoritmos.

Esse formato .hex é exclusivo do H2O e pode ser usado para inúmeros algoritmos dentro da plataforma, dado que ele já é otimizado para lidar com objetos esparsos e/ou colunas do tipo texto.

Para ver as estatísticas descritivas desse arquivo, basta usar o mesmo summary() do R.

# Let's see the summary
summary(airlines.hex)

Para o nosso experimento, vamos dividir a base de treino e teste na proporção 70%/30%.

Uma coisa necessária a se dizer nesse ponto é que devido ao fato do H2O ser uma plataforma projetada para Big Data é utilizado o método de amostragem probabilística. Isso se faz necessário (em termos computacionais), dado que em muitas vezes a operação de seleção/estratificação pode ser custoso.

# Construct test and train sets using sampling
# A small note is that H2O uses probabilistic splitting, witch means that resulting splits
# can deviate for the exact number. This is necessary when we're talking about a platform that 
# deals with big data. If you need a exact sampling, the best way is to do this in your RDBMS
airlines.split = h2o.splitFrame(data = airlines.hex,ratios = 0.70, seed = -1)

Após criar o objeto do tipo splitFrame, vamos alocar as partições para cada conjunto de dados, sendo que o objeto na primeira posição sempre será a nossa base de treinamento, e na segunda posição a nossa base de teste.

# Get the train dataframe(1st split object)
airlines.train = airlines.split[[1]]

# Get the test dataframe(2nd split object)
airlines.test = airlines.split[[2]]

Vamos sumarizar abaixo cada um desses frames para verificarmos a distribuição de voos cancelados:

# Display a summary using table-like in some sumarized way
h2o.table(airlines.train$Cancelled)
# Cancelled Count
# 1         0 29921
# 2         1   751

h2o.table(airlines.test$Cancelled)
# Cancelled Count
# 1         0 12971
# 2         1   335

Com as nossas amostras separadas, vamos agora escolher as variáveis que vão entrar no nosso modelo.

Primeiramente, vamos criar dois objetos para passar como parâmetro ao nosso algoritmo.

Como queremos prever se as partidas do voos estão atrasadas, então o objeto Y (variável dependente) será a variável IsDepDelayed (Se o voo de partida está atrasado) e o objeto X (variáveis independentes) serão todos os outros campos do conjunto de dados.

# Set dependent variable (Is departure delayed)
Y = "IsDepDelayed"
# Set independent variables
X = c("Origin", "Dest", "DayofMonth", "Year", "UniqueCarrier", "DayOfWeek", "Month", "DepTime", "ArrTime", "Distance")

Agora vamos realizar a criação do modelo usando GLM:

# Define the data for the model and display the results
airlines.glm <- h2o.glm(training_frame=airlines.train
                        ,x=X
                        ,y=Y
                        ,family = "binomial"
                        ,alpha = 0.5
                        ,max_iterations = 300
                        ,beta_epsilon = 0
                        ,lambda = 1e-05
                        ,lambda_search = FALSE
                        ,early_stopping = FALSE
                        ,nfolds = 0
                        ,seed = NULL
                        ,intercept = TRUE
                        ,gradient_epsilon = -1
                        ,remove_collinear_columns = FALSE
                        ,max_runtime_secs = 10000
                        ,missing_values_handling = c("Skip"))

O significado dos parâmetros do modelo são:

x: vetor que contém os nomes das variáveis independentes;

y: índice que contém a variável dependente;

training_frame: Um frame do H2O que contém das variáveis do modelo;

family: Especificação da distribuição do modelo que pode ser gaussiana, binomial, poisson, gamma, e tweedie. Uma ótima explicação de como esses parâmetros podem ser escolhidos está aqui nesse link;

alpha: Um número em [0, 1] especificando a mistura do parâmetro de regularização do elastic-net. Ele que dá o grau de mistura entre os regularizadores Lasso e Ridge. making alpha = 1 penalização via LASSO, alpha = 0 penalização via ridge;

max_iterations: Um inteiro não negativo que especifica o número máximo de interações do modelo;

beta_epsilon: Um inteiro não negativo que especifica a magnitude da diferença máxima entre as estimativas dos coeficientes através de sucessivas interações. Em outras palavras: É esse parâmetro que define a velocidade da convergência do modelo GLM;

lambda: Um parâmetro não negativo para encolhimento do valor da variável através da Elastic-Net, o qual multiplica P(α, β) na função objetivo. Quando lambda = 0, nenhuma penalização é aplicada e o modelo já fica ajustado;

lambda_search: Um valor lógico que indica se haverá algum critério de busca no espaço dos valores de lambda definidos por um parâmetro de mínimo e máximo;

early_stopping: Um valor que indica se haverá uma parada antecipada durante o lambda_search caso o fator de verosimilhança pare de ser alterado na medida que ocorram mais interações;

nfolds: Número de partições em Cross-Validation;

seed: Especifica a semente do random number generator (RNG) para Cross-Validation (garante a reprodutibilidade do experimento);

intercept: Termo constante do modelo que a grosso modo significa o grau de fatores endógenos do modelo;

gradient_epsilon: Critério de convergência. Converge se o gradiente da norma I-Infinito é abaixo de um determinado limite. Se lambda_search = FALSE e lambda = 0, o valor default do gradient_epsilon é igual a .000001, se não for, o valor default é .0001. Se lambda_search = TRUE, os valores condicionais acima são 1E-8 e 1E-6 respectivamente.

remove_collinear_columns: Se não houver nenhum tipo de fator de regularização aplicado, o modelo ignora colunas colineares (o coeficiente será 0);

max_runtime_secs: Número máximo permitido em segundos para o treinamento do modelo. Use 0 para desabilitar; e

missing_values_handling: Contra o que é feito com valores faltantes. Podem ser “MeanImputation” ou “Skip”. MeanImputation substituí os valores faltantes com a média para os atributos numéricos e categórico com a maior frequência. É aplicado durante o treinamento do modelo.

Notem que aqui o céu é o limite em temos de ajustes e/ou parametrizações. O ideal é ter o perfeito entendimento da mecânica de cada um dos parâmetros e utilizar a melhor combinação possível.

Com o nosso modelo ajustado, vamos ver algumas das estatísticas básicas desse modelo.

# View model information: training statistics, performance, important variables
summary(airlines.glm)

# Model Details:
#   ==============
#   
#   H2OBinomialModel: glm
# Model Key:  GLM_model_R_1484053333586_1 
# GLM Model: summary
# family  link                              regularization number_of_predictors_total number_of_active_predictors number_of_iterations  training_frame
# 1 binomial logit Elastic Net (alpha = 0.5, lambda = 1.0E-5 )                        283                         272                    5 RTMP_sid_a6c9_1
# 
# H2OBinomialMetrics: glm
# ** Reported on training data. **
#   
#   MSE:  0.2098326
# RMSE:  0.4580749
# LogLoss:  0.607572
# Mean Per-Class Error:  0.3720209
# AUC:  0.7316312
# Gini:  0.4632623
# R^2:  0.1602123
# Null Deviance:  41328.6
# Residual Deviance:  36240.45
# AIC:  36786.45
# 
# Confusion Matrix for F1-optimal threshold:
#   NO   YES    Error          Rate
# NO     5418  9146 0.627987   =9146/14564
# YES    1771 13489 0.116055   =1771/15260
# Totals 7189 22635 0.366047  =10917/29824
# 
# Maximum Metrics: Maximum metrics at their respective thresholds
# metric threshold    value idx
# 1                       max f1  0.363651 0.711915 294
# 2                       max f2  0.085680 0.840380 389
# 3                 max f0point5  0.539735 0.683924 196
# 4                 max accuracy  0.521518 0.673887 207
# 5                max precision  0.987571 1.000000   0
# 6                   max recall  0.040200 1.000000 398
# 7              max specificity  0.987571 1.000000   0
# 8             max absolute_mcc  0.521518 0.348709 207
# 9   max min_per_class_accuracy  0.513103 0.672412 212
# 10 max mean_per_class_accuracy  0.521518 0.674326 207

Aqui neste modelo já temos um resultado de 73,16% de AUC. Nada mal para um modelo que contém poucos ajustes.

Vamos analisar agora a importância de cada uma das variáveis no modelo:

# Get the variable importance of the models
h2o.varimp(airlines.glm)

# Standardized Coefficient Magnitudes: standardized coefficient magnitudes
# names coefficients sign
# 1 Origin.TLH     3.233673  NEG
# 2 Origin.CRP     2.998012  NEG
# 3 Origin.LIH     2.859198  NEG
# 4   Dest.LYH     2.766090  POS
# 5 Origin.KOA     2.461819  NEG
# 
# ---
#   names coefficients sign
# 278   Dest.JAN     0.000000  POS
# 279   Dest.LIT     0.000000  POS
# 280   Dest.SJU     0.000000  POS
# 281 Origin.LAN     0.000000  POS
# 282 Origin.SBN     0.000000  POS
# 283 Origin.SDF     0.000000  POS

Alguns valores de atributos são determinantes no atraso dos vôos de partida, principalmente se os aeroportos de origem são TLH (Tallahassee International Airport), CRP (Corpus Christi International Airport), LIH (Lihue Airport), e KOA (Kona International Airport).

Agora vamos usar a função predict para saber como o modelo realizou as classificações.

# Predict using GLM model
pred = h2o.predict(object = airlines.glm, newdata = airlines.test)

Após isso, vamos ver o resultado do nosso modelo.

# Look at summary of predictions: probability of TRUE class (p1)
summary(pred)

# predict   NO                YES              
# YES:9798  Min.   :0.01186   Min.   :0.02857  
# NO :3126  1st Qu.:0.33715   1st Qu.:0.37018  
# NA : 382  Median :0.48541   Median :0.51363  
#           Mean   :0.48780   Mean   :0.51220  
#           3rd Qu.:0.62886   3rd Qu.:0.66189  
#           Max.   :0.97143   Max.   :0.98814  
#           NA's   :382       NA's   :382 

Outra forma de usar o modelo GLM no H2O é realizar a escolha de parâmetros via Grid Search.

Grid Search nada mais é do que um método de otimização de parâmetros de um modelo de Machine Learning, sendo que em grande parte das vezes é feita uma combinação com inúmeros parâmetros e a combinação que obtiver o menor erro através de uma determinada função de erro.

O que vamos fazer agora é usar o GLM pra obter o melhor modelo de acordo com uma combinação de parâmetros específica.

Primeiramente, vamos construir uma lista com os parâmetros alpha (recapitulando, esse parâmetro faz uma combinação de Lasso e Ridge via Elastic-Net). Neste caso vamos passar uma lista que vai desde 0.0 até 1 com incremento de 0.05 em cada parâmetro.

# Construct a hyper-parameter space
alpha_opts = c(0,0.05,0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1)

Agora vamos criar uma lista com esses parâmetros para passar para a função de Grid posteriormente.

# List of hyperparameters
hyper_params_opt = list(alpha = alpha_opts)

Na função de Grid Search, vamos passar alguns parâmetros para o ajuste dos modelos como podemos ver abaixo.

# Grid object with hyperparameters
glm_grid <- h2o.grid("glm"
                     ,grid_id = "glm_grid_1"
                     ,x=X
                     ,y=Y
                     ,training_frame=airlines.train
                     ,hyper_params = hyper_params_opt
                     ,family = "binomial")

Essa etapa pode demorar bastante tempo dependendo do seu volume de dados, e do número de parâmetros escolhidos na lista de search.

Após a finalização do processamento, vamos ordernar a lista de modelos de acordo com o AUC.

# Sort grids by best performance (lower AUC). Little note: As we're dealing with classification
# in some probabilistc fashion, we'll use AUC as model selection metric.
# If the nature of the problem are cost sensitive (e.g. A delayed departure plane is much expensive for 
# the airport service than a delayed arrival) precision and recall can be the best choice
glm_sorted_grid <- h2o.getGrid(grid_id = "glm_grid_1", sort_by = "auc", decreasing = FALSE)

Para avaliar cada um dos modelos, podemos exibir a ordem dos modelos de acordo com o AUC.

#Print the models
print(glm_sorted_grid)

# H2O Grid Details
# ================
#   
#   Grid ID: glm_grid_1 
# Used hyper parameters: 
#   -  alpha 
# Number of models: 21 
# Number of failed models: 0 
# 
# Hyper-Parameter Search Summary: ordered by increasing auc
# alpha          model_ids                auc
# 1 [D@4800a43e glm_grid_1_model_1 0.7076911403181928
# 2 [D@66030470 glm_grid_1_model_2 0.7122987232329416
# 3 [D@6a4a43d3 glm_grid_1_model_3 0.7145455620514375
# 4 [D@17604a1a glm_grid_1_model_4  0.715989429818657
# 5 [D@21e1e99f glm_grid_1_model_5 0.7169797604977775
#                
# ---
# alpha           model_ids                auc
# 16 [D@78833412 glm_grid_1_model_16  0.720595118360825
# 17 [D@44d770f2 glm_grid_1_model_17 0.7207086912177467
# 18 [D@31669527 glm_grid_1_model_18 0.7208228330257134
# 19 [D@5b376f34 glm_grid_1_model_19 0.7209144533220885
# 20 [D@6acad45e glm_grid_1_model_20 0.7209885192412766
# 21 [D@237ad7de  glm_grid_1_model_0 0.7240682725570593

Com esses parâmetros, o melhor modelo é o glm_grid_1_model_0 que teve cerca de 72.40% de AUC. (Nota: Esse modelo está levemente pior do que o modelo padrão, dado que o conjunto de parâmetros do Grid está diferente do que o primeiro modelo).

 Para pegar o melhor modelo, basta executar o comando abaixo:

# Grab the model_id based in AUC
best_glm_model_id <- glm_grid@model_ids[[1]]
# The best model
best_glm <- h2o.getModel(best_glm_model_id)

Vejamos as características desse modelo:

# Summary
summary(best_glm)

# Model Details:
#   ==============
#   
#   H2OBinomialModel: glm
# Model Key:  glm_grid_1_model_0 
# GLM Model: summary
# family  link             regularization number_of_predictors_total number_of_active_predictors number_of_iterations  training_frame
# 1 binomial logit Ridge ( lambda = 7.29E-5 )                        283                         282                    3 RTMP_sid_a6c9_1
# 
# H2OBinomialMetrics: glm
# ** Reported on training data. **
#   
#   MSE:  0.2121424
# RMSE:  0.4605891
# LogLoss:  0.612699
# Mean Per-Class Error:  0.3833898
# AUC:  0.7240683
# Gini:  0.4481365
# R^2:  0.1494395
# Null Deviance:  42448.59
# Residual Deviance:  37585.41
# AIC:  38151.41
# 
# Confusion Matrix for F1-optimal threshold:
#   NO   YES    Error          Rate
# NO     4993  9601 0.657873   =9601/14594
# YES    1751 14327 0.108907   =1751/16078
# Totals 6744 23928 0.370110  =11352/30672
# 
# Maximum Metrics: Maximum metrics at their respective thresholds
# metric threshold    value idx
# 1                       max f1  0.373247 0.716243 296
# 2                       max f2  0.105583 0.846435 391
# 3                 max f0point5  0.551991 0.685249 194
# 4                 max accuracy  0.513313 0.665949 218
# 5                max precision  0.980714 1.000000   0
# 6                   max recall  0.048978 1.000000 399
# 7              max specificity  0.980714 1.000000   0
# 8             max absolute_mcc  0.548278 0.332916 196
# 9   max min_per_class_accuracy  0.524282 0.664324 211
# 10 max mean_per_class_accuracy  0.548278 0.666166 196
# 
# Gains/Lift Table: Extract with `h2o.gainsLift(&lt;model&gt;, &lt;data&gt;)` or `h2o.gainsLift(&lt;model&gt;, valid=&lt;T/F&gt;, xval=&lt;T/F&gt;)`
# 
# 
# 
# Scoring History: 
#   timestamp   duration iteration negative_log_likelihood objective
# 1 2017-01-10 11:11:07  0.000 sec         0             21224.29620   0.69198
# 2 2017-01-10 11:11:07  0.066 sec         1             18857.11178   0.61705
# 3 2017-01-10 11:11:07  0.094 sec         2             18795.11788   0.61562
# 4 2017-01-10 11:11:07  0.126 sec         3             18792.70362   0.61559
# 
# Variable Importances: (Extract with `h2o.varimp`) 
# =================================================
#   
#   Standardized Coefficient Magnitudes: standardized coefficient magnitudes
# names coefficients sign
# 1 Origin.MDW     1.915481  POS
# 2 Origin.HNL     1.709757  NEG
# 3 Origin.LIH     1.584259  NEG
# 4 Origin.HPN     1.476562  POS
# 5 Origin.AUS     1.439134  NEG
# 
# ---
#   names coefficients sign
# 278 Origin.PHX     0.009111  POS
# 279   Dest.PWM     0.008332  POS
# 280 Origin.GEG     0.008087  POS
# 281   Dest.BOS     0.005105  POS
# 282   Dest.MCI     0.003921  NEG
# 283   Dest.CHA     0.000000  POS

Para realizar previsões com esse modelo, basta apenas instanciar esse novo objeto e usar a função predict como está abaixo:

# Get model and put inside a object
model = best_glm

# Prediction using the best model
pred2 = h2o.predict(object = model, newdata = airlines.test)

# Summary of the best model
summary(pred2)

# predict    NO                YES              
# YES:10368  Min.   :0.01708   Min.   :0.05032  
# NO : 2938  1st Qu.:0.33510   1st Qu.:0.39258  
#            Median :0.47126   Median :0.52781  
#            Mean   :0.47526   Mean   :0.52474  
#            3rd Qu.:0.60648   3rd Qu.:0.66397  
#            Max.   :0.94968   Max.   :0.98292  

Se após isso, você quiser desligar o cluster do H2O basta usar o comando shutdown.

# Shutdown the cluster 
h2o.shutdown()

# Are you sure you want to shutdown the H2O instance running at http://localhost:54321/ (Y/N)? Y
# [1] TRUE

Com isso finalizamos esse post/tutorial de como usar o R com o H2O.ai.

Ao longo das próximas semanas, vamos trazer alguns tutoriais e destrinchar um pouco o poder desse H2O.

Forte abraço!

 

 

 

 

 

 

 

 

 

 

Anúncios
Implementação de GLM com Grid Search no R usando o H2O.ai como backend

2 comentários sobre “Implementação de GLM com Grid Search no R usando o H2O.ai como backend

Deixe o seu comentário inteligente e educado! :o)

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s