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

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:

[sourcecode language=”r”] #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! ;) [/sourcecode]

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

[sourcecode language=”r”] # 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”) } [/sourcecode] Em seguida vamos fazer o download e instalação de todos os pacotes dos quais o H2O tem alguma dependência direta ou indireta. [sourcecode language=”r”] # 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”) } [/sourcecode] Em seguida faremos a instalação da lib do H2O e o instanciamento da lib no R Studio. [sourcecode language=”r”] # 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”))) [/sourcecode] [sourcecode language=”r”] # Load library library(h2o) [/sourcecode]

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.

[sourcecode language=”r”] # 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) [/sourcecode] Agora vamos ver as informações do nosso cluster: [sourcecode language=”r”] # 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) [/sourcecode]

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.

[sourcecode language=”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” [/sourcecode] [sourcecode language=”r”] # 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”) [/sourcecode]

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.

[sourcecode language=”r”] # Let’s see the summary summary(airlines.hex) [/sourcecode]

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.

[sourcecode language=”r”] # 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) [/sourcecode]

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.

[sourcecode language=”r”] # Get the train dataframe(1st split object) airlines.train = airlines.split[[1]] # Get the test dataframe(2nd split object) airlines.test = airlines.split[[2]] [/sourcecode]

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

[sourcecode language=”r”] # 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 [/sourcecode]

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.

[sourcecode language=”r”] # Set dependent variable (Is departure delayed) Y = “IsDepDelayed” [/sourcecode] [sourcecode language=”r”] # Set independent variables X = c(“Origin”, “Dest”, “DayofMonth”, “Year”, “UniqueCarrier”, “DayOfWeek”, “Month”, “DepTime”, “ArrTime”, “Distance”) [/sourcecode] Agora vamos realizar a criação do modelo usando GLM: [sourcecode language=”r”] # 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”)) [/sourcecode]

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.

[sourcecode language=”r”] # 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 [/sourcecode]

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:

[sourcecode language=”r”] # 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 [/sourcecode]

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.

[sourcecode language=”r”] # Predict using GLM model pred = h2o.predict(object = airlines.glm, newdata = airlines.test) [/sourcecode] Após isso, vamos ver o resultado do nosso modelo. [sourcecode language=”r”] # 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 [/sourcecode]

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.

[sourcecode language=”r”] # 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) [/sourcecode] Agora vamos criar uma lista com esses parâmetros para passar para a função de Grid posteriormente. [sourcecode language=”r”] # List of hyperparameters hyper_params_opt = list(alpha = alpha_opts) [/sourcecode] Na função de Grid Search, vamos passar alguns parâmetros para o ajuste dos modelos como podemos ver abaixo. [sourcecode language=”r”] # 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”) [/sourcecode]

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.

[sourcecode language=”r”] # 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) [/sourcecode]

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

[sourcecode language=”r”] #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 [/sourcecode]

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:

[sourcecode language=”r”] # Grab the model_id based in AUC best_glm_model_id <- glm_grid@model_ids[[1]] [/sourcecode] [sourcecode language=”r”] # The best model best_glm <- h2o.getModel(best_glm_model_id) [/sourcecode] Vejamos as características desse modelo: [sourcecode language=”r”] # 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(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)` # # # # 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 [/sourcecode] Para realizar previsões com esse modelo, basta apenas instanciar esse novo objeto e usar a função predict como está abaixo: [sourcecode language=”r”] # 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 [/sourcecode] Se após isso, você quiser desligar o cluster do H2O basta usar o comando shutdown. [sourcecode language=”r”] # 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 [/sourcecode]

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!