Análise de Múltipla Correspondência no R para o problema de Churn

Via Data Science Plus

Analytical challenges in multivariate data analysis and predictive modeling include identifying redundant and irrelevant variables. A recommended analytics approach is to first address the redundancy; which can be achieved by identifying groups of variables that are as correlated as possible among themselves and as uncorrelated as possible with other variable groups in the same data set. On the other hand, relevancy is about potential predictor variables and involves understanding the relationship between the target variable and input variables.
Multiple correspondence analysis (MCA) is a multivariate data analysis and data mining tool for finding and constructing a low-dimensional visual representation of variable associations among groups of categorical variables. Variable clustering as a tool for identifying redundancy is often applied to get a first impression of variable associations and multivariate data structure.
The motivations of this post are to illustrate the applications of: 1) preparing input variables for analysis and predictive modeling, 2) MCA as a multivariate exploratory data analysis and categorical data mining tool for business insights of customer churn data, and 3) variable clustering of categorical variables for the identification of redundant variables.

Análise de Múltipla Correspondência no R para o problema de Churn

Deep Dive com Gradient Boosting Machine com H2O + R (Mais Grid Search!)

Dando sequência a alguns tutoriais sobre o uso do R como linguagem de programação junto H2O como backend de processamento e memória (duas principais limitações do R) vamos falar um pouco de Gradient Boosting Machine e usar uma base de dados de crédito de um banco fictício chamado “Layman Brothers”.

Gradient Boosting Machine é um meta-algoritmo de aprendizado supervisionado que é geralmente utilizado em problemas de classificação e regressão. O principio algorítmico por trás do GBM é a produção de previsões/classificações derivadas de modelos preditivos fracos (Weak Learners), em especial árvores de decisão essas que por sua vez combinadas via ensemble learning para redução de vieses dos algoritmos.

Essas previsões são geradas através da combinação da meta-heurística de gradiente descendente para otimização paramétrica face a minimização de uma função de custo (loss function), e do Boosting que é combinação de diversos classificadores fracos (Weak Learners) em série para (ou meta-classificador) para combinação de resultados desses algoritmos.

Como podemos supor, com essa combinação heurística de algoritmos, em especial dos weak learners (que dão uma robustez substancial ao modelo) é de se esperar uma determinada insensibilidade á distribuição de cauda longa que pode ser espessa e detonar as suas previsões (e.g. distribuição da renda mundial em que poucos (20%) tem muito dinheiro e muitos (80%) tem pouco) , outliers (i.e. eventos extremos, também conhecidos como cisnes negros), além de uma boa resposta a não-linearidade. (Nota: Se você não entendeu nada do que está aqui, uma boa pedida são dois livros do Nassim Taleb que são Black Swan (A lógica do cisne negro) e Antifragile (Antifrágil)).

Como dito anteriormente, a base de dados que será usada aqui é de um banco fictício chamado “Layman Brothers”, que é uma alusão simpática ao Lehman Brothers; e o nosso objetivo é ter um sistema de crédito um pouco mais confiável do que o deles o que não é uma tarefa que demande muita inteligência ou stamina intelectual. (Nota: Essa base é originalmente do repositório do UCI, mas estou rebatizando para dar um tom cênico mais descontraído aqui no post).

A nossa base de dados de créditos tem as seguintes colunas:

  • ID: Número da transação
  • LIMIT_BAL: Crédito concedido em dólares
  • SEX: Sexo (1 = masculino; 2 = feminino).
  • EDUCATION: Nível escolar d@ cliente (1 = ensino médio; 2 = universidade; 3 = ensino superior completo; 4 = outros)
  • MARRIAGE: Estado civil (1 = casad@; 2 = solteir@; 3 = outros).
  • AGE: Idade d@ cliente
  • PAY_X: Histórico do pagamento passado. Foi rastreado o pagamento passado mensal (de abril até setembro de 2005) da seguinte forma: PAY_1 o status de repagamento do mês de setembro de 2005, PAY_2: o status do repagamento mês de agosto de 2005, etc. A escala de medida do repagamento é :-1 = Pago em dia, 1 = pago com um mês de atraso, 2 = pagamento atrasado por 2 meses, 8 = pagamento atrasado por 8 meses, etc.
  • BILL_AMTX: Montante do saldo ainda não amortizado dos meses anteriores. BILL_AMT1 = Saldo ainda não amortizado em setembro de 2005, BILL_AMT2 = saldo ainda não amortizado em agosto de 2005, etc.
  • PAY_AMTX: Montante pago anteriormente (em dólares) relativos ao mês anterior. PAY_AMT1 = valor pago em setembro de 2005, PAY_AMT2 = valor pago em agosto de 2005, etc.
  • DEFAULT: Se @ cliente deixou de pagar o empréstimo no mês seguinte.

Base de dados apresentada, vamos ao código.

Primeiramente, se você não instalou o H2O via R ou está com a versão desatualizada, é só executar esse código abaixo que ele vai remover a versão antiga, instalar todas as dependências, e instalar o H2O:

# 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") }

# 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") }

# 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")))

Agora vamos carregar a biblioteca e iniciar o nosso cluster (que nesse caso ainda estará no meu notebook) com o tamanho máximo de memória de 8 gigas, e vai usar todos os processadores (-1):

# Load library
library(h2o)


# Start instance with all cores
h2o.init(nthreads = -1, max_mem_size = "8G")

# Info about cluster
h2o.clusterInfo()

# Production Cluster (Not applicable because we're using in the same machine)
#localH2O <- h2o.init(ip = '10.112.81.210', port =54321, nthreads=-1) # Server 1
#localH2O <- h2o.init(ip = '10.112.80.74', port =54321, nthreads=-1) # Server 2

Cluster iniciado, vamos buscar os nossos dados que estão no repositório remoto do Github e na sequência vamos carregar no nosso objeto .hex (extensão do H2O):

# URL with data
LaymanBrothersURL = "https://raw.githubusercontent.com/fclesio/learning-space/master/Datasets/02%20-%20Classification/default_credit_card.csv"

# Load data 
creditcard.hex = h2o.importFile(path = LaymanBrothersURL, destination_frame = "creditcard.hex")

Com os dados carregados, vamos realizar a transformação das variáveis categóricas, e em seguida vamos ver o sumário dessas variáveis:

# Convert DEFAULT, SEX, EDUCATION, MARRIAGE variables to categorical
creditcard.hex[,25] <- as.factor(creditcard.hex[,25]) # DEFAULT
creditcard.hex[,3] <- as.factor(creditcard.hex[,3]) # SEX
creditcard.hex[,4] <- as.factor(creditcard.hex[,4]) # EDUCATION
creditcard.hex[,5] <- as.factor(creditcard.hex[,5]) # MARRIAGE

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

Como podemos ver pelo summary() temos algumas estatísticas descritivas básicas interessantes sobre essa base de dados, como:

screen-shot-2017-01-15-at-12-17-16-pm

  • A maioria dos empréstimos foram feitos por pessoas que se declararam do sexo feminino (60%);
  • 63% de todos os empréstimos foram feitos para a população classificada como universitária ou que tem curso superior completo;
  • Há um equilíbrio entre o estado civil em relação aos empréstimos concedidos;
  • Com um terceiro quartil de 41 e uma média e medianas bem próximas (35 e 34), podemos ver que grande parte dos empréstimos foram feitos por pessoas na idade adulta que estão na meia idade; e
  • Temos muitas pessoas que pegaram empréstimos altos (acima de 239 mil dólares), porém, a média do valor concedido é de 167 mil dólares.

Óbvio que caberiam mais algumas análises de perfil, correlações, e até mesmo alguns gráficos para exemplificar melhor a composição demográfica dessa base, mas como esse não é o objetivo desse post, fica aberto para que algum dos 5 leitores desse site blog faça isso e compartilhe.

Com essas análises feitas, vamos dividir a nossa base nos conjuntos de treinamento, teste e validação usando o comando splitFrame:

# We'll get 3 dataframes Train (60%), Test (20%) and Validation (20%)
creditcard.split = h2o.splitFrame(data = creditcard.hex
                                  ,ratios = c(0.6,0.2)
                                  ,destination_frames = c("creditcard.train.hex", "creditcard.test.hex", "creditcard.validation.hex")
                                  ,seed = 12345)


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

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

# Get the validation dataframe(3rd split object)
creditcard.validation = creditcard.split[[3]]

Para checarmos a real proporção de cada base, podemos usar o comando table para ver a composição de cada base de dados (e principalmente ver se elas estão balanceadas):

# See datatables from each dataframe
h2o.table(creditcard.train$DEFAULT)

# DEFAULT Count
# 1       0 14047
# 2       1  4030

h2o.table(creditcard.test$DEFAULT)

# DEFAULT Count
# 1       0  4697
# 2       1  1285

h2o.table(creditcard.validation$DEFAULT)

# DEFAULT Count
# 1       0  4620
# 2       1  1321

Agora vamos criar dois objetos para passar ao nosso algoritmo: um objeto para definir quem será a nossa variável dependente (Y) e outro para definir as nossas variáveis independentes (X):

# Set dependent variable
Y = "DEFAULT"

# Set independent variables
X = c("LIMIT_BAL","EDUCATION","MARRIAGE","AGE"
      ,"PAY_0","PAY_2","PAY_3","PAY_4","PAY_5","PAY_6"
      ,"BILL_AMT1","BILL_AMT2","BILL_AMT3","BILL_AMT4","BILL_AMT5","BILL_AMT6"
      ,"PAY_AMT1","PAY_AMT3","PAY_AMT4","PAY_AMT5","PAY_AMT6")

# I intentionally removed sex variable from the model, to avoid put any gender bias inside the model. Ethics first guys! 😉

Os mais atentos podem verificar que eu removi a variável SEX. Fiz isso intencionalmente dado que não vamos colocar nenhum tipo de viés discriminatório no modelo (Atenção amigos: esse é um bom tempo para considerar seriamente essas questões de discriminação/ética em modelos de Machine Learning como etnia, gênero, etc).

Agora com esses objetos prontos, vamos treinar o nosso modelo:

# Train model
creditcard.gbm <- h2o.gbm(y = Y
                          ,x = X
                          ,training_frame = creditcard.train
                          ,validation_frame = creditcard.validation                      
                          ,ntrees = 100
                          ,seed = 12345
                          ,max_depth = 100
                          ,min_rows = 10
                          ,learn_rate = 0.2
                          ,distribution= "bernoulli"
                          ,model_id = 'gbm_layman_brothers_model'
                          ,build_tree_one_node = TRUE
                          ,balance_classes = TRUE
                          ,score_each_iteration = TRUE
                          ,ignore_const_cols = TRUE
                          )

Explicando alguns desses parâmetros:

  • x: Vetor que contém os nomes das variáveis independentes do modelo;
  • y: índice ou objeto que representa a variável dependente do modelo;
  • training frame: Um objeto de dados do H2O (H2OFrame) que contém as variáveis do modelo;
  • validation frame: Um objeto de dados do H2O (H2OFrame) que contém as variáveis do modelo para validação do modelo. Se estiver vazia os dados de treinamento são usados por padrão;
  • ntrees: Um inteiro não negativo que define o número de árvores. O valor default é 50;
  • seed: Semente dos números aleatórios a serem gerados. É usado para reprodutibilidade amostral;
  • max depth: Valor definido pelo usuário do número máximo da profundidade das árvores. O valor default é 5;
  • min rows: O número mínimo de linhas a serem designadas para cada nó terminal. O padrão é 10;
  • learn rate: Um inteiro que define a taxa de aprendizado do modelo. Vai de 0.1 até 1.0;
  • distribution: Escolhe uma distribuição de probabilidade entre AUTO, bernoulli, multinomial, gaussian, poisson, gamma ou tweedie. O default é AUTO;
  • model id: ID único que identifica o modelo. Se não especificado é gerado automaticamente;
  • build tree one node: Especifica se o modelo será processado em um nó apenas. Isso serve para evitar overhead de rede e com isso menos CPUs são usadas no processo. É ideal para pequenos datasets, e o default é FALSE;
  • balance classes: Faz o balanceamento de classes do conjunto de treinamento, caso os dados estejam com subamostragem ou desbalanceados. O default é falso;
  • score each iteration: Um binário que indica se haverá o processo de scoring durante cada interação do modelo. O default é falso; e
  • ignore const cols: Um binário que indica se colunas com constantes serão ignoradas. O Default é TRUE.

Alguns conselhos práticos de quem já sofreu (muito) na pele para parametrizar GBM que você não vai ter do seu professor na faculdade:

a) O H2O oferece e a opção validation_frame, porém, se você for mais purista o ideal é checar na etapa de prediction e ver o bias do modelo através da análise dos erros (sim gente, vai ter que rolar estatística aqui, ok?). Isso além de dar um ajuste mais fino, te dá o maior entendimento dos erros modelo. Se fosse em minas, o pessoal lá diria que isso faz bem pra saúde e forma o caráter. Faça o mesmo.;

b) Tenha bastante parcimônia para ajustar o número ideal de árvores (ntrees) dado que isso eleva demais o custo computacional (processamento + memória) do modelo. Via de regra, eu gosto de usar intervalos de 50 árvores para cada step até o limite de 300; e assim que eu chego em um meio termo eu vou ajustando na unha via grid search até chegar em uma árvore que eu tenha um bom desempenho sem overfitting. Isso é necessário pois grande parte das vezes você tem uma elevação ridícula de até 8 horas no tempo de treinamento pra ganhar no máximo 0.01 no AUC, ou uma redução de 0.005% nos falsos positivos. Em resumo: Vai com calma no ajuste. Faz bem pra saúde e forma o caráter; e além do mais economiza mais de 20 dólares na Amazon pra treinar um modelo caso você esteja usando máquinas on-demand fora da sua infra;

c) É o seed que vai garantir que os seus números estão corretos quando você for passar para alguém fazer o code review ou mesmo antes do deployment. Então use sempre que puder por questões óbvias de reprodutibilidade;

d) O parâmetro max depth costuma ser o que eu chamo de cemitério do malandro em Machine Learning. Isso devido ao fato de que qualquer iniciante em seu primeiro contato com esse parâmetro vai colocar o maior número possível em geral quase o mesmo número de instâncias da base de treinamento (isso é quando o malandro não coloca cross-validation pra coisa ficar ainda mais bonita) o que deixa a árvore mais específica e leva na maioria das vezes aquele overfittingTem iniciantes que conseguem a proeza de fazer overfitting mesmo usando max depth com leave-one-out cross validation. (Pequena dica empírica: pessoalmente eu nunca consegui resultados bacanas com uma profundidade de níveis que excedam 0.005% do número de registros no conjunto de treinamento (100/((30000/100)*70 =0.005%). Ainda estou tentando saber se isso está correto ou não, mas ao menos pra mim funciona bem;

e)  Quanto menor o valor do min rows, mais específica será a árvore e pode ocorrer que ela generalize menos. Por isso muita parcimônia com esse parâmetro;

f) Desnecessário dizer que um número muito pequeno pode influenciar no tempo de processamento e convergência do modelo, e um número alto pode cair em um mínimo local e estragar todo o trabalho. Dica prática: tá com pouco tempo? Vai de 0.35 até 0.75 com incremento de 0.1. Tá com tempo de sobra? Vai de 0.1 até 0.5 com incremento de 0.03;

g) Realmente vale a pena gastar um pouco de neurônios para entender melhor as distribuições de probabilidade (distribution) para escolher a correta. Se você não tiver tempo, escolha a AUTO e seja feliz;

h) A não ser que você esteja enfrentando uma situação de concorrência de rede e de processamento, o parâmetro build tree one node sempre deve estar desligado;

i) Se você está usando o parâmetro balance classes significa que o seu trabalho de amostragem está um lixo e você precisa da ferramenta pra fazer algo básico pode não ser o mais correto. Eu recomendo fortemente uma seriedade no processo de amostragem que é o coração de qualquer treinamento de machine learning. Caso sejam situações amostrais muito esquisitas (e.g. modelagem de sistemas de combate á fraudes, classificador de reclamações em Call Center, et cetera) ou por falta de tempo, vale a pena usar esse parâmetro (Dica prática: caso haja uma situação de desbalanceamento muito grave de classes (algo na proporção 9:1) o ideal é esquecer as outras métricas de avaliação de modelos e ir direto para o coeficiente de matthews que é bem mais consistente para lidar com esse tipo de caso);

j) Se você está usando o parâmetro ignore const cols é porque o seu trabalho de pré-processamento (Feature Extraction e Feature Engineering) está um lixo pode não estar sendo o melhor.

Modelo treinado e parâmetros explicados, vamos ver a performance do modelo usando os dados de validação:

# See algo performance
h2o.performance(creditcard.gbm, newdata = creditcard.validation)

# H2OBinomialMetrics: gbm

# MSE:  0.1648487
# RMSE:  0.4060157
# LogLoss:  0.8160863
# Mean Per-Class Error:  0.3155595
# AUC:  0.7484422
# Gini:  0.4968843

# Confusion Matrix for F1-optimal threshold:
#   0    1    Error        Rate
# 0      3988  632 0.136797   =632/4620
# 1       653  668 0.494322   =653/1321
# Totals 4641 1300 0.216294  =1285/5941

# We have an AUC of 74,84%, not so bad!

Com esse modelo tivemos um AUC de 74,84%. Razoável, considerando que usamos um conjunto de parametrizações simples.

A seguir, vamos conferir a importância de cada uma de nossas variáveis:

# Variable importance
imp <- h2o.varimp(creditcard.gbm)

head(imp, 20)

# Variable Importances: 
#   variable relative_importance scaled_importance percentage
# 1  EDUCATION        17617.437500          1.000000   0.380798
# 2   MARRIAGE         9897.513672          0.561802   0.213933
# 3      PAY_0         3634.417480          0.206297   0.078557
# 4        AGE         2100.291992          0.119217   0.045397
# 5  LIMIT_BAL         1852.831787          0.105170   0.040049
# 6  BILL_AMT1         1236.516602          0.070187   0.026727
# 7   PAY_AMT5         1018.286499          0.057800   0.022010
# 8  BILL_AMT3          984.673889          0.055892   0.021284
# 9  BILL_AMT2          860.909119          0.048867   0.018608
# 10  PAY_AMT6          856.006531          0.048589   0.018502
# 11  PAY_AMT1          828.846252          0.047047   0.017915
# 12 BILL_AMT6          823.107605          0.046721   0.017791
# 13 BILL_AMT4          809.641785          0.045957   0.017500
# 14  PAY_AMT4          771.504272          0.043792   0.016676
# 15  PAY_AMT3          746.101196          0.042350   0.016127
# 16 BILL_AMT5          723.759521          0.041082   0.015644
# 17     PAY_3          457.857758          0.025989   0.009897
# 18     PAY_5          298.554657          0.016947   0.006453
# 19     PAY_4          268.133453          0.015220   0.005796
# 20     PAY_2          249.107925          0.014140   0.005384

Nesse modelo podemos ver que o nível educacional tem um papel essencial na definição de quem vai entrar em default (38%), seguindo do estado civil (21%) e fechando com o pagamento anterior relativo ao mês de setembro de 2008 (7%) e da idade do tomador de crédito e o saldo emprestado (4%).

Em outras palavras: essas variáveis acima respondem por 74% do comportamento de crédito.

Com isso algumas questões hipóteses (Hx) e ações (Ax) podem ser tomadas pelo Layman Brothers:

H1: O nível educacional está muito relacionado com o default,  isso acontece de forma positiva ou não em relação à inadimplência?

H2: Será que universitários que tradicionalmente são pessoas com menos poder aquisitivo tem maiores dificuldades (ou facilidades) para o pagamento?

H3: De que forma o estado civil influencia na capacidade de pagamento do crédito emprestado?

H4: Porque o saldo não amortizado exerce efeito tão grande em relação às outras variáveis financeiras?

H5: Porque a pontualidade no pagamento não é tão determinante, com exceção da primeira parcela?

H6: O perfil educacional influencia o quanto em relação à capacidade de pagamento?

A1: De acordo com a escolaridade, ter diferentes taxas de juros para empréstimos.

A2: Ter ações de cobrança efetivas/intensas já no primeiro mês de atraso.

A3: Ter linhas de crédito mais específicas para cada perfil educacional com taxas e saldos correspondentes ao risco de default.

A4: Entender e criar linhas de financiamento de acordo com cada objetivo de acordo com o estado civil (e.g. entender se o gasto é para investimento (voltado para a geração de mais receita como cursos, maquinário, ou outros fatores que aumentem a produtividade; ou para despesas como consumo, contas de inúmeras naturezas, outros empréstimos, et cetera) .

Adiante, podemos agora usar o nosso modelo treinado para fazer previsões:

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

# See predictions
head(pred, 5)

# predict        p0           p1
# 1       0 0.9990856 0.0009144487
# 2       0 0.9945627 0.0054373206
# 3       0 0.9997726 0.0002273775
# 4       0 0.9968271 0.0031728833
# 5       0 0.9991758 0.0008242144

Agora, vamos para um ajuste mais fino no nosso modelo com o objetivo de melhorar o nosso AUC (que é atualmente de 74,84%), e para isso vamos usar Grid Search.

Primeiramente vamos gerar uma lista de valores para os nossos hiper-parâmetros (hyper parameters) do modelo GBM. Os parâmetros que vamos usar serão ntrees (número de árvores), max_depth (profundidade das árvores) e learn_rate (taxa de aprendizado). Após isso vamos jogar dentro de uma meta lista que vamos usar para ajustar o nosso objeto de grid.

# Set hyparameters (Did not work using sequence. :o( )
ntrees_list <- list(50,100,150,200)

max_depth_list <- list(1,2,3,4,5,6,7,8,9,10)

learnrate_list <- list(.10,.20,.30,.40,.50,.60,.70,.80,.90)
# Full list of hyper parameters that will be used
hyper_parameters <- list(ntrees = ntrees_list
                         ,max_depth = max_depth_list
                         ,learn_rate = learnrate_list)

# See hyparameters lists
hyper_parameters

Ou seja, teremos uma combinação com 50, 100, 150 e 200 árvores, níveis de profundidade da árvore indo de 1 até 10 e taxa de aprendizado indo de 0.10 até 0.90.

Uma pequena experiência da trincheira deste escriba que não foi muito inteligente é ter uma boa combinação de números de parâmetros na meta lista em relação com a capacidade de processamento disponível para fazer o treinamento.

Isso se faz necessário pois como abaixo vamos usar a estratégia cartesiana para o nosso critério de busca (i.e. vamos usar todas as combinações paramétricas possíveis) vamos ter o seguinte cenário:

ntrees = 4
max_depth = 10
learn_rate = 9

Logo teremos 4 * 10 * 9 = 360 modelos/combinações!

Ou seja: Pode levar bastante tempo para processar (no meu caso levou 11m34min pra acabar, e houve uma porção de erros do H2O por incapacidade de processamento).

Após o processamento do grid vamos ordenar os modelos do melhor para o pior usando o AUC:

# sort the grid models by decreasing AUC
sortedGrid <- h2o.getGrid("depth_grid", sort_by="auc", decreasing = TRUE)    
# Let's see our models
sortedGrid

# H2O Grid Details
# ================
  
# Grid ID: depth_grid 
# Used hyper parameters: 
# -  learn_rate 
# -  max_depth 
# -  ntrees 
# Number of models: 380 
# Number of failed models: 2940 

# Hyper-Parameter Search Summary: ordered by decreasing auc
# learn_rate max_depth ntrees            model_ids                auc
# 1        0.1         6    100 depth_grid_model_200 0.7811807105334736
# 2        0.1         6     50   depth_grid_model_5 0.7811440893197138
# 3        0.2         3    150 depth_grid_model_264 0.7809025695475355
# 4        0.2         3    100 depth_grid_model_174  0.780834324645831
# 5        0.1         6    200 depth_grid_model_380 0.7808292451933633

Agora, vamos pegar o melhor modelo (com menor AUC) e vamos ver algumas das suas características:

# Summary
summary(best_glm)

# Model Details:
# ==============
  
# H2OBinomialModel: gbm
# Model Key:  depth_grid_model_200 
# Model Summary: 
#   number_of_trees number_of_internal_trees model_size_in_bytes min_depth max_depth mean_depth
# 1             100                      100               52783         6         6    6.00000
# min_leaves max_leaves mean_leaves
# 1         12         56    36.93000

# H2OBinomialMetrics: gbm
# ** Reported on training data. **
  
# MSE:  0.1189855
# RMSE:  0.3449427
# LogLoss:  0.3860698
# Mean Per-Class Error:  0.2593832
# AUC:  0.8371354
# Gini:  0.6742709

# Confusion Matrix for F1-optimal threshold:
# 0    1    Error         Rate
# 0      12424 1623 0.115541  =1623/14047
# 1       1625 2405 0.403226   =1625/4030
# Totals 14049 4028 0.179676  =3248/18077

Esse nosso modelo tem 100 árvores, uma profundidade de 6 níveis, e em média 37 instâncias em cada nó folha.

Como podemos ver tivemos um AUC de 83,71%, ou 11% de melhoria em comparação com o antigo AUC que foi de 74,84% em menos de 12 minutos.

Um fato curioso é que olhando a importância das variáveis novamente com esse modelo temos os seguintes resultados:

# Variable importance (again...)
imp2 <- h2o.varimp(best_glm)

head(imp2, 20)

# Variable Importances: 
#   variable relative_importance scaled_importance percentage
# 1      PAY_0         2040.270508          1.000000   0.358878
# 2      PAY_2          902.637390          0.442411   0.158772
# 3  LIMIT_BAL          385.425659          0.188909   0.067795
# 4        AGE          274.609589          0.134595   0.048303
# 5  BILL_AMT1          209.715469          0.102788   0.036888
# 6      PAY_3          168.518372          0.082596   0.029642
# 7  EDUCATION          150.365280          0.073699   0.026449
# 8  BILL_AMT2          146.754837          0.071929   0.025814
# 9      PAY_5          139.303482          0.068277   0.024503
# 10  PAY_AMT5          139.206543          0.068229   0.024486
# 11 BILL_AMT5          133.963348          0.065660   0.023564
# 12     PAY_4          124.926552          0.061230   0.021974
# 13  PAY_AMT6          123.267151          0.060417   0.021682
# 14 BILL_AMT6          114.012253          0.055881   0.020054
# 15  PAY_AMT1          112.402290          0.055092   0.019771
# 16     PAY_6          108.483795          0.053171   0.019082
# 17 BILL_AMT3          103.207893          0.050585   0.018154
# 18  PAY_AMT3           97.335411          0.047707   0.017121
# 19 BILL_AMT4           90.403320          0.044309   0.015902
# 20  MARRIAGE           61.917801          0.030348   0.010891

Ou seja, se antigamente o nível educacional e o estado civil tiveram uma participação importante, nesse modelo (com melhor AUC) a pontualidade, o montante de crédito concedido e a idade exercem mais influência.

Com esse melhor modelo, podemos fazer as nossas previsões e salvar em um arquivo .csv para upload em algum sistema ou isso pode ser feito via API via requisição.

# Get model and put inside a object
model = best_glm

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

# Frame with predictions
dataset_pred = as.data.frame(pred2)

# Write a csv file
write.csv(dataset_pred, file = "predictions.csv", row.names=TRUE)

 E após finalizado todo o trabalho, podemos desligar o nosso cluster:

# 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

Bem pessoal como vocês podem ver, usar um modelo usando Gradient Boosting Machine no R não é nenhum bicho de 7 cabeças no H2O, basta um pouquinho de parcimônia na parametrização que tudo dá certo.

Se tiverem dúvidas deixem o seu comentário inteligente e educado aqui nos comentários ou me mandem por e-mail.

Forte abraço!

 

Deep Dive com Gradient Boosting Machine com H2O + R (Mais Grid Search!)

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!

 

 

 

 

 

 

 

 

 

 

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

Principais soluções do H2O.ai

Agora que já sabemos um pouco sobre essa solução, vamos entender um pouco do ecossistema de soluções do H2O, e ver as principais características e aplicações de cada uma.

H2O

Essa plataforma é o carro chefe da empresa, o qual eles apostam tanto na versão Desktop para aplicação de Machine Learning quando também na versão para processamento distribuído para altos volumes de dados.

Essa versão tem alguns algoritmos prontos on the shelf como boosting, regressão linear e logística, algoritmos baseados em árvores, e alguns algoritmos que utilizam gradiente como método de otimização. Nada muito complexo, mas bem funcional.

Essa versão é ideal para usar se você quer conhecer um pouco mais da ferramenta e não quer gastar muito tempo instalando ou configurando coisas antes de sair aplicando os algoritmos, ou mesmo para um teste inicial das funções de processamento distribuído em cluster. 

Abaixo, um pouco da arquitetura da solução:

h2oarch

Sparkling Water

Essa solução tem como principal característica utilizar os próprios algoritmos, mas com a vantagem de usar todas as features de processamento distribuído e integração do Spark. Nesta solução todas as tarefas de computação também podem ser feitas dentro do Spark usando Scala e com uma interface via Web.

Essa solução é a mais recomendada para construção de aplicações de Machine Learning seja para microserviços ou até mesmo para embutir dentro de uma plataforma/sistema toda a parte algorítmica e computacional.

h2ospark h2ospark2

Deep Water

O Deep Water é a solução voltada para implementação de Deep Learning usando otimização computacional com GPUs com frameworks como o Tensor Flow, Theano, Caffe entre outros.

Neste caso, a plataforma do H2O será a interface onde serão incorporados todos os parâmetros de treinamento do modelo (cross validation, amostragem, critério de parada, hiperparametrização, etc) e o backend com o Tensor Flow, Theano etc. faz o processamento utilizando GPUs.

Steam

O Steam é uma plataforma que realiza todo o link entre os modelos de machine learning usando o H2O e também propriedades de desenvolvimento para incorporar modelos de Machine Learning em aplicações, tudo isso de forma colaborativa, algo muito similar ao Domino.

A principal vantagem do Steam é que ele abstraí toda a parte de engenharia por trás da tarefa de incorporar modelos e machine learning em produção como infraestrutura, auto-scale de infra estrutura de acordo com a carga de requisições, bem com algumas tarefas de Data Science como retreinamento de modelos; além de reduzir e muito os custos/investimentos de TI.

steam

Agora que sabemos quais são os principais produtos do H2O, em breve teremos alguns posts com alguns tutoriais explorando um pouco mais essa ferramenta.

Links úteis

Documentação Técnica

Principais soluções do H2O.ai

Porque você REALMENTE deveria considerar o H2O.ai como uma das ferramentas do seu stack de Machine Learning?

Quando falamos de ferramentas de machine learning logo vem a cabeça a tríade Tensor Flow,  Scikit-Learn e Spark MLLib.
Contudo, há uma ferramenta que vem discretamente ganhando espaço que é o H2O.ai.
Essa ferramenta nasceu originalmente em 2011 em que o seu time teve como principal objetivo democratizar e tornar escalável machine learning através de uma plataforma mais visual e que tivesse uma boa experiência para os usuários, independente do seu nível técnico.
Algumas características do H2O.ai:
Abaixo alguns vídeos sobre o H2O.ai em ação:
H2O.ai para detecção de fraudes
Customer Churn usando H2O.ai

Em alguns posts futuros falaremos um pouco sobre questões de arquitetura, e mergulharemos em alguns tutoriais.
Porque você REALMENTE deveria considerar o H2O.ai como uma das ferramentas do seu stack de Machine Learning?

Previsão de Séries Temporais usando XGBoost – Pacote forecastxgb

Para quem já teve a oportunidade de trabalhar com previsão de variáveis categóricas em Machine Learning sabe que o XGBoost é um dos melhores pacotes do mercado, sendo largamente utilizado em inúmeras competições no Kaggle.

A grande diferença feita pelo Peter Ellis foi realizar algumas adaptações para incorporar algumas variáveis independentes através do parâmetro xreg ao modelo preditivo de séries temporais.

Para quem trabalha com análise de séries temporais, esse trabalho é muito importante até porque o forecastxgb  tríade Média-Móvel/ARIMA (ARMA)/(S)ARIMA em que tanto estatísticos/Data Miners/Data Scientists ficam presos por comodidade ou falta de meios.

Um exemplo da utilização do pacote está abaixo:

# Install devtools to install packages that aren't in CRAN
install.packages("devtools")

# Installing package from github 
devtools::install_github("ellisp/forecastxgb-r-package/pkg") 

# Load the libary
library(forecastxgb)

# Time Series Example
gas

# Model
model <- xgbts(gas)

# Summary of the model
summary(model)

# Forecasting 12 periods 
fc <- forecast(model, h = 12)

# Plot
plot(fc)
Previsão de Séries Temporais usando XGBoost – Pacote forecastxgb

Regras de Associação no Combate ao Crime no Brasil – Por Flávio Clesio

Não é segredo que vivemos em um país que está lamentavelmente listado como um dos mais violentos do mundo, em que temos 21 cidades entre as 50 mais violentas do mundo.

Se levarmos em consideração que em alguns estados do Brasil a taxa de esclarecimento de homicídios é de menos de 15%, temos um problema grave na segurança pública de maneira geral.

Mas se o crime por si só é a externalidade/consequência de diversos fatores como aspectos socioeconômicos, educacionais, e de ambiente (para citar alguns); o policiamento preventivo, isto é, o ataque direto das causas deve ser o norte de qualquer política de segurança pública; e é sobre isso que se trata esse post.

Através de uma base de dados sintética (extremamente desbalanceada) de dados que poderiam se extraídos de qualquer distrito policial ou mesmo de algum departamento de segurança pública vamos fazer a aplicação de uma técnica de Machine Learning chamada Association Rules ou Regras de Associação para formulação de políticas de segurança pública.

Como vamos iniciar a nossa análise sem nenhum tipo de hipótese e o objetivo é ver o conhecimento intrínseco da base de dados, vamos fazer um aprendizado não-supervisionado.

Contexto colocado, vamos para o R Studio.

Antes de mais nada, vamos limpar o nosso console e na sequência importar a biblioteca _arules_ no R.

# Clear the console
cat("\014") 

# Load arules library
library(arules)

Com a nossa lib importada, vamos realizar carregar a nossa base de dados no R.

# Load Dataset - This is a synthetic dataset, some inconsistencies can be found... as in the real world! 😉
crimes <- read.csv("https://raw.githubusercontent.com/fclesio/learning-space/master/Datasets/01%20-%20Association%20Rules/Crimes.csv")

Com a nossa base de dados importada, vamos agora ver as colunas da base de dados:

# See variables types
str(crimes)

Temos uma base de dados com 780 crimes com 16 variáveis que são:

  • Crime_id: Número da ocorrência.
  • Zona: Zona da cidade que a ocorrência foi registrada.
  • Periculosidade: Nível de periculosidade da área de acordo com mapeamentos anteriores.
  • Viaturas_Apoio_Forca_Tatica: Indica se tem reforços da Força Tática disponíveis na região. A grosso modo, Força Tática é um destacamento que ‘flutua’ na circulação e tem como objetivo atender a ocorrência mais próxima que a primeira força de atendimento não conseguir solucionar.
  • Patrulhamento: Indica se há patrulhamento pelas viaturas do batalhão mais próximo.
  • Policiamento_Ostensivo: Informa se há operações acontecendo com mais frequência. Um maior nível de policiamento aumenta também o enfrentamento policial, o que pode acarretar em mais óbitos de ambas as partes.
  • Apoio_GCM: Esse campo indica se a Guarda Civil Metropolitana circula na área da ocorrência. Diferente da PM ou da Civil; a GCM apesar das mesmas prerrogativas policiais tem uma função mais voltada à segurança patrimonial de alguns pontos da cidade e escolas.
  • Arma_Fogo: Houve arma de fogo na ocorrência.
  • Qtde_Vitimas: Número de vitimas da ocorrência.
  • Possui_DP: Tem Distrito Policial nas redondezas.
  • Tipo_Policiamento: O tipo de policiamento local influencia também na prática de crimes, em que áreas menos assistidas pelos batalhões podem ter um maior índice de violência em relação à áreas que não tem esse tipo de atendimento com um maior volume de viaturas e policiais.
  • Area_Residencial: Indica se é área residencial ou não.
  • Ocorrencia_Atendida_15_Minutos: Sabe-se que quanto mais demorado for o atendimento da ocorrência, as chances de resolução do evento criminoso em questão converge para a não solução efetiva. Esse campo indica se a ocorrência foi prontamente atendida.
  • Pericia: Informa se houve a perícia para auxiliar nas investigações.
  • Possui_UPP: Mostra se tem Unidade de Polícia Pacificadora instalada nas redondezas do evento criminoso
  • Iluminacao: Informa se o local do crime tinha iluminação ou não.

Com um pouco mais de trabalho na parte de data wrangling é possível criar mais algumas variáveis, mas como o foco aqui será de fazer uma primeira análise básica vamos seguir com algumas estatísticas básicas em relação a cada uma das variáveis da base de dados.

# See basic stats
summary(crimes)

Uma olhada rápida no sumário podemos tirar as seguintes informações:

  • Em relação ao tipo de policiamento vemos que temos um patrulhamento feito por diversos batalhões na metade dos casos, o que indica que há uma dispersão de batalhões (i.e. nenhum batalhão cuida específicamente de uma área) nessas áreas mais críticas;
  • 69% das localidades tem algum tipo de Distrito Policial nas redondezas do evento criminoso;
  • Apenas 1% dos locais não tem nenhum tipo de viatura da Força Tática;
  • Quase 70% das ocorrências acontecem com 3 ou mais pessoas. Isso nos faz a crer que pode ter dois tipos de vieses que são a) os crimes coletivos (e.g. arrastões, chacinas, etc.) são os mais reportados ou estão acontecendo de fato, ou b) os crimes menores como roubo de celular não estão sendo notificados;
  • Mais de metade das áreas dos eventos aconteceram em áreas consideradas com o nível de periculosidade altíssimo;
  • Pelos dados a Zona Sul está em chamas;
  • Com 88% dos locais com iluminação, pode ser que essa variável não seja mais a responsável pela inibição dos eventos criminosos;
  • Com apenas 49% dos casos periciados, as atividades de investigação podem em algum momento ter problemas de robustez de provas; e
  • Sendo o índice de patrulhamento em 90% e o policiamento ostensivo 83% nas regiões criminosas há indícios de que os criminosos estão indiferentes à presença policial.

São insights interessantes que nos auxiliam em algumas hipóteses como as levantadas. Neste ponto a análise descritiva básica pode ajudar e muito no início de uma política de segurança efetiva, e não pode ser descartada em nenhum momento. No nosso caso vamos avançar para a parte aplicada de machine learning que é o foco desse post.

Vamos criar um objeto utilizando o algoritmo a priori escolhendo um suporte de 55% (isso é, os conjuntos de itens estarão no mínimo em 55% de todas as transações) e com a confiança de 80% (isso é, a grosso modo é como dizer que há a probabilidade de 80% das transações encontradas serão derivadas dos itens contidos na base do suporte).

Não há uma regra única em relação à utilização do suporte e confiança, mas trazendo aqui uma heurística que eu sempre uso é que se o problema tem um volume de transações muito grande, o suporte será a métrica que vai dar mais robustez para descartar muitas regras redundantes, já se o problema tem uma caracterização com um menor volume de instâncias e se há uma necessidade de entrar mais no detalhe das regras e elas precisam de um grau de assertividade naturalmente maior (e.g. uma campanha de cartas em que há um custo operacional envolvido) aí utiliza-se a confiança. Mas lembrando, tudo isso depende do problema em questão.

Nesse caso, como os recursos a serem empregados são caros e escassos (viaturas e pessoal militar), então partiremos para uma confiança de 80%.

# The most important thing it's get a good support and a great confidence (This is a wishful thinking, in real world EVERYTHING can be happen!)
# The main task here is get a set a most meaningful way, so, shity rules will be discarted.
crimes_rulez <- apriori(crimes, parameter = list(minlen=4,maxlen=10, supp=0.55, conf=0.80))

Com o objeto criado, vamos agora reduzir para duas casas decimais as métricas de avaliação dessas regras para fins de simplicidade.

# Reduce the number of digits in measures of quality of the model
quality(crimes_rulez) <- round(quality(crimes_rulez), digits=2)

Apenas para fins de visualização, vamos criar alguns objetos de regras ordenando pela confiança e suporte para que possamos identificar de maneira mais efetiva as regras mais interessantes e acionáveis.

# Sort by Confidence and Support and Lift
model_srt_confidence <- sort(crimes_rulez, by="confidence")

model_srt_support <- sort(crimes_rulez, by="support")

Vamos agora inspecionar as nossas regras, de acordo com a confiança e o suporte:

# Check models (#66 rules)
inspect(model_srt_confidence)

inspect(model_srt_support)

Com esse comando acima, temos algumas regras:

      
   lhs                                     rhs                                  support confidence lift

1  {Apoio_GCM=Nao,                                                                                     
    Area_Residencial=Sim,                                                                              
    Iluminacao=Sim}                     => {Patrulhamento=Sim}                     0.57       0.98 1.09

63 {Area_Residencial=Sim,                                                                              
    Ocorrencia_Atendida_15_Minutos=Sim,                                                                
    Iluminacao=Sim}                     => {Zona=Zona_Sul}                         0.55       0.82 0.96

49 {Patrulhamento=Sim,                                                                                 
    Area_Residencial=Sim,                                                                              
    Ocorrencia_Atendida_15_Minutos=Sim,                                                                
    Iluminacao=Sim}                     => {Policiamento_Ostensivo=Sim}            0.58       0.89 1.08
  • A regra 1 indica que o patrulhamento está funcionando nos locais em que a PM não tem o apoio da GCM em áreas iluminadas de bairros residenciais;
  • A regra de número 63 mostra um primeiro padrão em relação a Zona Sul, em que com 82% de confiança essas ocorrências estão ligadas a lugares com malha residencial, iluminadas e a ocorrência foi feita dentro dos primeiros 15 minutos;
  • Já a regra 49 que o policiamento ostensivo está ocorrendo com 89% de confiança em áreas já patrulhadas em uma malha residencial iluminada em que o evento foi atendido em 15 minutos ou menos.

Com essas poucas regras acima podemos observar que o trabalho de primeiro atendimento está sendo feito dado que grande parte dessas ocorrências (dado o nível de suporte e confiança) estão sendo atendidas em um tempo de 15 minutos.

Outro ponto a ser considerado é que nas áreas residenciais estão ocorrendo muitos eventos criminosos. Isso pode ser um indicativo de que nas áreas de comércio e indústrias podem estar com baixa atratividade criminosa ou mesmo o policiamento nessas áreas está mais efetivo.  Uma estratégia aqui seria ver esses dois aspectos, e caso a segunda hipótese estiver correta replicar as ações de combate ao crime.

Já a Zona Sul que concentra grande parte das ocorrências, parece estar em chamas mesmo com um alto nível de patrulhamento, iluminação e policiamento ostensivo. Neste caso há a hipótese de que são o maior patrulhamento e um policiamento mais ostensivo além de ter uma ação mais rápida – que converge como o fato do atendimento de ocorrências demorar menos de 15 minutos –  ela tem o poder de jogar luz em uma demanda reprimida de registros policiais (i.e. crimes que não seriam registrados por ausência do poder público, ou mesmo por descaso do estado nunca entrariam nas estatísticas oficiais).

Em um segundo caso pode ser mesmo que uma facção criminosa esteja disputando controle direto com a polícia através de confrontos armados, ou mesmo atentados e roubos em áreas que são vigiadas pela polícia; o que sugere uma estratégia de inteligência para desarticular de dentro para fora milícias, grupos criminosos ou quadrilhas free lance que estejam levando terror a essas áreas.

Como são muitas regras e o objetivo é encontrar o máximo de regras específicas, vamos minerar algumas dessas regras usando campos específicos. Sendo assim vamos explorar um pouco mais das ocorrências que apresentaram as características de terem ocorrido na Zona Sul, em lugares iluminados, casos que as ocorrências foram atendidas em menos de 15 minutos ou que tiveram policiamento ostensivo.

# See specific rules of zones (I already mined some rules, but feel free to discover)
Zona_Sul <- subset(crimes_rulez, subset = rhs %pin% "Zona=Zona_Sul")

Iluminacao_Sim <- subset(crimes_rulez, subset = rhs %pin% "Iluminacao=Sim")

Ocorrencia_Atendida_15_Minutos_Sim <- subset(crimes_rulez, subset = rhs %pin% "Ocorrencia_Atendida_15_Minutos=Sim")

Policiamento_Ostensivo <- subset(crimes_rulez, subset = rhs %pin% "Policiamento_Ostensivo=Sim")

Esses objetos criados, tem somente as regras em que o RHS (right-hand-side) foram esses campos descritos no parágrafo anterior.

Vamos inspecionar esses objetos.

# Specific Rules
inspect(Zona_Sul)

inspect(Iluminacao_Sim)

inspect(Ocorrencia_Atendida_15_Minutos_Sim)

inspect(Policiamento_Ostensivo)

Seguindo a mesma lógica da análise acima, podemos ter os seguintes conjuntos de regras específicas:

Zona_Sul
8  {Patrulhamento=Sim,                                                            
    Policiamento_Ostensivo=Sim,                                                   
    Iluminacao=Sim}                     => {Zona=Zona_Sul}    0.66       0.88 1.03

9  {Patrulhamento=Sim,                                                            
    Area_Residencial=Sim,                                                         
    Iluminacao=Sim}                     => {Zona=Zona_Sul}    0.66       0.85 0.99


Iluminacao_Sim
   lhs                                     rhs              support confidence lift
12 {Patrulhamento=Sim,                                                             
    Policiamento_Ostensivo=Sim,                                                    
    Area_Residencial=Sim,                                                          
    Ocorrencia_Atendida_15_Minutos=Sim} => {Iluminacao=Sim}    0.58       0.96 1.09

13 {Zona=Zona_Sul,                                                                 
    Patrulhamento=Sim,                                                             
    Policiamento_Ostensivo=Sim,                                                    
    Area_Residencial=Sim}               => {Iluminacao=Sim}    0.61       0.95 1.08    


Ocorrencia_Atendida_15_Minutos_Sim
  lhs                             rhs                                  support confidence lift
5 {Zona=Zona_Sul,                                                                             
   Area_Residencial=Sim,                                                                      
   Iluminacao=Sim}             => {Ocorrencia_Atendida_15_Minutos=Sim}    0.55       0.81 1.08

7 {Patrulhamento=Sim,                                                                         
   Policiamento_Ostensivo=Sim,                                                                
   Area_Residencial=Sim,                                                                      
   Iluminacao=Sim}             => {Ocorrencia_Atendida_15_Minutos=Sim}    0.58       0.83 1.11


Policiamento_Ostensivo
   lhs                                     rhs                          support confidence lift
10 {Patrulhamento=Sim,                                                                         
    Area_Residencial=Sim,                                                                      
    Ocorrencia_Atendida_15_Minutos=Sim,                                                        
    Iluminacao=Sim}                     => {Policiamento_Ostensivo=Sim}    0.58       0.89 1.08

11 {Zona=Zona_Sul,                                                                             
    Patrulhamento=Sim,                                                                         
    Area_Residencial=Sim,                                                                      
    Iluminacao=Sim}                     => {Policiamento_Ostensivo=Sim}    0.61       0.93 1.12       

No caso dessas regras, tudo indica que mesmo com as ocorrências sendo atendidas em menos de 15 minutos em áreas residenciais iluminadas, a Zona Sul está com um problema grave de segurança.

A estratégia nesse caso seria investir de forma intensa em atividades de inteligência, dado que pelo que vimos pelos dados há uma situação de combate intenso com policia na rua e crimes acontecendo mesmo em situações em que os bandidos supostamente estariam intimidados.

Conclusão
Usando uma base sintética, vimos que é possível com algumas linhas de código e estatística descritiva básica ter insights sobre uma base de dados mesmo sem um conhecimento prévio.

No caso da nossa análise, vimos que pode haver uma situação de confronto clara em áreas residenciais mesmo com patrulhamento da polícia e policiamento ostensivo acontecendo.

Isso leva a crer que esses crimes estejam sendo reportados pelos aspectos acima, ou mesmo por uma maior disposição dos criminosos em agir em áreas controladas pela polícia. Isso pode indicar que o trabalho da polícia nesse caso não está sendo eficiente, e que uma estratégia de investigação através de atividades de inteligência devem ser colocadas em campo.

Regras de Associação no Combate ao Crime no Brasil – Por Flávio Clesio

Novo pacote do R – forecastHybrid

Direto do Peter Stats Blog

The new forecastHybrid package for R by David Shaub and myself provides convenient access to ensemble time series forecasting methods, in any combination of up to five of the modelling approaches from Hyndman’s forecast package. These are auto-selected autoregressive integrated moving average; exponential smoothing state space (both ets and tbats); feed forward neural network with a single hidden layer and lagged inputs; and forecasts of loess-based seasonal decomposition.

Background and motivation
In an earlier post I explored ways that might improve on standard methods for prediction intervals from univariate time series forecasting. One of the tools I used was a convenience function to combine forecasts from Rob Hyndman’s ets and auto.arima functions. David Shaub (with a small contribution from myself) has now built and published an R package forecastHybrid that expands on this idea to create ensembles from other forecasting methods from Hyndman’s forecast package.

The motivation is to make it easy to improve forecasts, both their point estimates and their prediction intervals. It has been well known for many years that taking the average of rival forecast methods improves the performance of forecasts. This new R package aims to make it as easy for people to do this as to fit the individual models in the first place.

Novo pacote do R – forecastHybrid

Gradient Boosted Trees Model no Microsoft R Server

Direto do Revolutions Blog

R is an open source, statistical programming language with millions of users in its community. However, a well-known weakness of R is that it is both single threaded and memory bound, which limits its ability to process big data. With Microsoft R Server (MRS), the enterprise grade distribution of R for advanced analytics, users can continue to work in their preferred R environment with following benefits: the ability to scale to data of any size, potential speed increases of up to one hundred times faster than open source R.

In this article, we give a walk-through on how to build a gradient boosted tree using MRS. We use a simple fraud data data set having approximately 1 million records and 9 columns. The last column “fraudRisk” is the tag: 0 stands for non-fraud and 1 stands for fraud. The following is a snapshot of the data.

Gradient Boosted Trees Model no Microsoft R Server

Previsão do número de homicídios no Rio de Janeiro (Sim, eu voltei a postar!!!)

Um pequeno disclaimer: Depois de muito relutar, estou voltando com as atividades do blog. Isso significa de maneira prática, que apesar do clipping ser o grande driver dos acessos, sinto que há uma lacuna em relação a posts técnicos aqui na blogosfera brasileira.

Dessa forma, vou fazer mais posts técnicos com assuntos ligados à Machine Learning, Data Mining e afins.

Disclaimer feito, bora pro post original.

*******

Fala pessoal, tudo bem?

Um dos principais problemas em utilizar processadores numéricos como o Excel ou mesmo o LibreOffice é que se perde um pouco a flexibilidade de parametrizar os modelos de previsão.

Além disso, há um problema grave nesses processadores de reprodutibilidade, dado que pode haver inúmeras versões, planilhas voando pra todo o lado via e-mail, e por aí vai.

Nesse post eu vou utilizar dados reais da quantidade de homicídios que houveram no Rio de Janeiro de 1991 até 2011.

Os dados foram extraídos e consolidados pelo CESeC – Centro de Estudos de Segurança e Cidadania e o link está no final do script.

Antes de mais nada, vamos carregar o package Forecast:

# Load library
if(!require(forecast)){
 install.packages("forecast")
 library(forecast)
}

Depois disso, vamos realizar a carga dos dados. Nesse caso eu fiz uma transformação na planilha original e passei os dados para csv.

# Dataset
rj_murder &lt;- read.csv("https://raw.githubusercontent.com/fclesio/learning-space/master/Datasets/rj-homicidios-1991-2011.csv")

Vamos dar uma olhada nos dados para nos certificar que tudo está OK.

# Check data
head(rj_murder)

Com os dados OK, para incluir um conjunto de dados no package de Forecast e ARIMA, temos que antes de mais nada, pegar todos os datapoints e construir uma série através de um objeto de timeseries como no bloco de código abaixo:

# Transform in time series dataset all metrics inside the dataset 
rj_murder_total &lt;- ts(rj_murder$total, start=c(1991, 1), end=c(2011, 12), frequency=12)
rj_murder_capital &lt;- ts(rj_murder$capital, start=c(1991, 1), end=c(2011, 12), frequency=12)
rj_murder_baixada &lt;- ts(rj_murder$baixada, start=c(1991, 1), end=c(2011, 12), frequency=12)
rj_murder_interior &lt;- ts(rj_murder$interior, start=c(1991, 1), end=c(2011, 12), frequency=12)

Para entendermos o que aconteceu, vamos destrinchar erssa função TS:

ts() função para a construção de série temporal

rj_murder$total conjunto de datapoints do dataframe da coluna total

start=c(1991, 1) início da série é janeiro de 1991

end=c(2011, 12) o fim da série é dezebro de 2011

frequency=12 como estamos falando de dados mensais, a frequencia da série é de 12 datapoints por ano

Agora que entendemos o que fizemos, vamos dar uma olhada nos nossos objetos ts() e ver se os mesmos estão corretos.

#Plot series
par(mfrow=c(2,2))

plot(rj_murder_total
 ,main = "Total de Homicidios no RJ de 1991-2011"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 )


plot(rj_murder_capital
 ,main = "Total de Homicidios no RJ de 1991-2011 (Capital)"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
)


plot(rj_murder_baixada
 ,main = "Total de Homicidios no RJ de 1991-2011 (Baixada)"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
)


plot(rj_murder_interior
 ,main = "Total de Homicidios no RJ de 1991-2011 (Interior)"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
)

 

Homicidios RJ

 

Agora que vimos que as nossas séries estão corretas, vamos criar alguns objetos usando a função auto.arima().

# Fit ARIMA models
fit_rj_murder_total &lt;- auto.arima(rj_murder_total)
fit_rj_murder_capital &lt;- auto.arima(rj_murder_capital) 
fit_rj_murder_baixada &lt;- auto.arima(rj_murder_baixada)
fit_rj_murder_interior &lt;- auto.arima(rj_murder_interior)

Essa função retorna o melhor modelo ARIMA de acordo com os valores de AIC, AICc e BIC.

Mas deixando um pouco a teoria de lado, vamos ver como ficaram os ajustes dos modelos ARIMA de acordo com cada uma das séries.

#Plot ARIMA Models
par(mfrow=c(2,2))
plot(forecast(fit_rj_murder_total,h=12)
 ,main = "Total de Homicidios no RJ de 1991-2011 \n Previsão usando ARIMA"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5)

plot(forecast(fit_rj_murder_capital,h=12)
 ,main = "Total de Homicidios no RJ de 91-2011 (Capital) \n Previsão usando ARIMA"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5)

plot(forecast(fit_rj_murder_baixada,h=12)
 ,main = "Total de Homicidios no RJ de 91-2011 (Baixada) \n Previsão usando ARIMA"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5)

plot(forecast(fit_rj_murder_interior,h=12)
 ,main = "Total de Homicidios no RJ de 91-2011 (Interior) \n Previsão usando ARIMA"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5)

 

ARIMA

 

Um bom resultado, mas agora vamos ver o que conseguimos fazer usando o package forecast default, com o horizonte de previsão de 12 meses:

#Forecasting using confidence intervals in 80%, 95% and 99% with 12 months ahead
ahead_rj_murder_total &lt;- forecast(rj_murder_total, level=c(80, 95, 99), h=12)
ahead_rj_murder_capital &lt;- forecast(rj_murder_capital, level=c(80, 95, 99), h=12)
ahead_rj_murder_interior &lt;- forecast(rj_murder_interior, level=c(80, 95, 99), h=12)
ahead_rj_murder_baixada &lt;- forecast(rj_murder_baixada, level=c(80, 95, 99), h=12)

Modelos ajustados, vamos plotar os resultados:

#Plot forecasting
par(mfrow=c(2,2))
plot(ahead_rj_murder_total
 ,main = "Total de Homicidios no RJ de 91-2011 (Total) \n Previsão usando Forecast package"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 ,shadecols="oldstyle")

plot(ahead_rj_murder_capital
 ,main = "Total de Homicidios no RJ de 91-2011 (Capital) \n Previsão usando Forecast package"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 ,shadecols="oldstyle")

plot(ahead_rj_murder_baixada
 ,main = "Total de Homicidios no RJ de 91-2011 (Baixada) \n Previsão usando Forecast package"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 ,shadecols="oldstyle")

plot(ahead_rj_murder_interior
 ,main = "Total de Homicidios no RJ de 91-2011 (Interior) \n Previsão usando Forecast package"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 ,shadecols="oldstyle")

Forecast Package Predictions

Podemos ver que tanto usando o modelo ARIMA quanto com o uso do package forecast vimos que haverá uma queda no número de homicídios. Com isso podem ser elaboradas (usando mais dados) políticas públicas de segurança.

Refs:

[1] – http://www.inside-r.org/packages/cran/forecast/docs/auto.arima

[2] – http://www.ucamcesec.com.br/wordpress/wp-content/uploads/2011/04/RJ_2011_evolucao.xls

[3] – http://www.statmethods.net/advstats/timeseries.html

 

Previsão do número de homicídios no Rio de Janeiro (Sim, eu voltei a postar!!!)

Análise Exploratoria de Dados utilizando gráficos

É fato que a inconsistência de dados acaba com qualquer tipo de modelagem em Data Mining.

Dessa forma, ANTES de qualquer experimento com data mining é sempre desejável que se faça uma análise exploratória de dados utilizando estatísticas descritivas, gráficos, formulação de hipóteses para uma definição clara de quais técnicas serão utilizadas.

Esse post do Eli Bressert apresenta um tutorial muito bom usando R para a utilização de gráficos na análise exploratória de dados.

Análise Exploratoria de Dados utilizando gráficos

Comparação entre R e Python utilizando Florestas Aleatórias e Classificação

Neste post do blog do Yhat tem o código, os dados e os resultados.

Pessoalmente gosto muito da abordagem dos autores em comparação de classificadores usando as métricas de Acurácia, Erro Quadrático Médio e  para regressão e tempo de treinamento.

Para projetos curtos de avaliação de uma série de classificadores essas medidas são suficientes para dar uma linha de base. Essas medidas podem auxiliar na escolha de quais modelos estão com melhor convergência e podem indicar um melhor tratamento dos dados em termos de quais variáveis são pertinentes ao modelo escolhido.

 

 

Comparação entre R e Python utilizando Florestas Aleatórias e Classificação