Introdução

Na primeira parte desta série, começamos analisando um banco de dados de churn customer com o pacote tidymodels. Nele discutimos a criação de etapas de pré-processamento dos dados com o pacote recipes, a criação de bases de treinamento e teste com rsample, a estimação de modelos com o parsnipe a criação de workflows com workflow. Por fim estimamos um modelo de random forest com alguns hiperparâmetros fixos que obteve um F-1 de 0,85.

Agora vamos utilizar cross-validation para determinar os valores de hiperparâmetros que minimizam algum critério de performance do modelo. Para isto, vamos utilizar os pacotes tune e dials do tidymodels.

Tuning e Cross-validation: tune e dials

Vamos começar importando os dados e declarando a mesma receita utilizada na parte 1 desta série.

library(tidyverse)
library(tidymodels)
library(skimr)
library(knitr)
library(doFuture)

data(wa_churn, package = "modeldata")

set.seed(123)

# base de treinamento e teste
tbl_treinamento_teste <- initial_split(data = wa_churn, prop = 0.8)

receita_simples <- recipe(churn ~ ., data = training(tbl_treinamento_teste)) %>% 
  step_impute_linear(all_numeric()) %>% 
  step_dummy(all_nominal(), -all_outcomes()) %>% 
  step_normalize(all_numeric())

Agora, em vez de utilizar um valor padrão para os hiperparâmetros do modelo, podemos ajustar diferentes valores utilizando a função tune() do pacote tune. Assim, deixamos que o procedimento de cross-validation decida o melhor valor dos hiperparâmetros do nosso modelo.

Para encontrar os valores ótimos dos hiperparâmetros podemos utilizar otimização bayesiana (função tune_bayes()) quando o problema de tuning for complexo. Contudo, para problemas mais simples como este, um grid search é suficiente (grid_regular). Assim, vamos construir uma combinação de valores dos hiperparâmetros, de modo que em cada fold do cross-validation serão estimados múltiplos modelos com diferentes hiperparâmetros.

Antes disto, precisamos construir nossa base de cross-validation. O pacote rsample permite construir um objetivo de cross-validation com n folds. Para minimizar o tempo de execução dos modelos, vamos construir apenas 4 folds diferentes. Em cada fold temos uma base de análise (equivalente ao treinamento) e de assessment (equivalente ao teste).

A função vfold_cv também aceita um argumento strata, que faz com que a amostragem aleatória seja conduzida de modo estratificado. No nosso exemplo, faremos uma amostragem aleatória simples. É possível observar que em todos os folds a proporção de churns é semelhante.

cv_folds <- vfold_cv(wa_churn, v = 4)

map_dbl(cv_folds$splits, ~ mean(as.data.frame(.x)$churn == "Yes"))

## [1] 0.2673230 0.2667550 0.2635365 0.2638652

Agora precisamos definir novamente nosso modelo, mas agora sem valores fixos para mtry, trees e min_n. Na parte 1, haviamos definido o modelo como:

modelo_rf <-
  rand_forest(mtry = 3, trees = 200, min_n = 30) %>% 
  set_mode("classification") %>% 
  set_engine("ranger")
modelo_rf

## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = 3
##   trees = 200
##   min_n = 30
## 
## Computational engine: ranger

Agora, vamos sinalizar com o uso de tune() que diferentes valores dos hiperparâmetros devem ser utilizados.

modelo_rf_tuning <-
  rand_forest(mtry = tune(), trees = tune(), min_n = tune()) %>% 
  set_mode("classification") %>% 
  set_engine("ranger")

modelo_rf_tuning

## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = tune()
##   min_n = tune()
## 
## Computational engine: ranger

E criar um workflow diferente, agora com um modelo que aceita hiperparâmetros não fixos.

workflow_rf_tuning <- workflow() %>% 
  add_model(modelo_rf_tuning) %>% 
  add_recipe(receita_simples)

workflow_rf_tuning

## == Workflow ====================================================================
## Preprocessor: Recipe
## Model: rand_forest()
## 
## -- Preprocessor ----------------------------------------------------------------
## 3 Recipe Steps
## 
## * step_impute_linear()
## * step_dummy()
## * step_normalize()
## 
## -- Model -----------------------------------------------------------------------
## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = tune()
##   min_n = tune()
## 
## Computational engine: ranger

Para utilizar a função grid_regular, precisamos informar os valores possíveis de cada hiperparâmetro. Para o randomForest, trees e min_n tem valores previamente definidos para o intervalo [1,2000] e [2,40], respectivamente. Contudo, mtry não possui um intervalo definido por padrão.

workflow_rf_tuning %>% parameters %>% {.$object}

## [[1]]
## # Randomly Selected Predictors (quantitative)
## Range: [1, ?]
## 
## [[2]]
## # Trees (quantitative)
## Range: [1, 2000]
## 
## [[3]]
## Minimal Node Size (quantitative)
## Range: [2, 40]

Para tanto, vamos atualizar nosso modelo com os valores máximos e mínimos dos hiperparâmetros utilizando update.parameters do pacote dials. Para mtry vamos definir seus valores mínimos e maximos como 1 e 30. trees é um hiperparâmetro quantitativo sem relação com a base de dados, seus valores máximos e mínimos já são definidos por padrão como 1 e 2000. O mesmo ocorre para min_n, que possui valor entre 2 e 40. Vamos guardar estas informações no objeto rf_parametros.

rf_parametros <- workflow_rf_tuning %>% 
  parameters() %>% 
  update(mtry = mtry(range = c(1L, 30L)),
         trees = trees(),
         min_n = min_n())

rf_parametros$object

## [[1]]
## # Randomly Selected Predictors (quantitative)
## Range: [1, 30]
## 
## [[2]]
## # Trees (quantitative)
## Range: [1, 2000]
## 
## [[3]]
## Minimal Node Size (quantitative)
## Range: [2, 40]

Os valores do grid poderiam ser definidos manualmente com o uso de uma função como expand.grid:

expand.grid(mtry = c(1,30), 
            trees = c(1,2000),
            min_2 = c(2,40)) 

##   mtry trees min_2
## 1    1     1     2
## 2   30     1     2
## 3    1  2000     2
## 4   30  2000     2
## 5    1     1    40
## 6   30     1    40
## 7    1  2000    40
## 8   30  2000    40

Mas grid_regular torna esta tarefa mais conveniente. Assim, passamos o argumento levels = 2 que indica que serão utilizados 2 valores igualmente esparçados dos três hiperparâmetros. Poderiamos ainda produzir valores aleatórios do grid utilizando a função random_grid().

rf_grid <- grid_regular(rf_parametros, levels = 2)
rf_grid

## # A tibble: 8 x 3
##    mtry trees min_n
##   <int> <int> <int>
## 1     1     1     2
## 2    30     1     2
## 3     1  2000     2
## 4    30  2000     2
## 5     1     1    40
## 6    30     1    40
## 7     1  2000    40
## 8    30  2000    40

A tabela acima mostra que temos \(\text{level}^m = 2 \times 2 \times 2\) combinações de hiperparâmetros. Isso significa que a realização de cross-validation e do tuning de hiperparâmetros aumenta enormente o número de modelos que serão ajustados. Temos 4 folds e 8 combinações possíveis de hiperparâmetros, totalizando 32 modelos.

Para aumentar a velocidade com que o processo de estimação ocorre, podemos estimar os modelos em paralelo com o uso do pacote doFuture. Para tanto, será necessário a instalações de alguns pacotes adicionais como Rmpi e snow, que dependem da instalação do Microsoft MPI que pode ser encontrado aqui. Caso ocorra dificuldades na instalação dos pacotes acima, é possível ignorar os códigos abaixo. Neste caso, os modelos podem ser ajustados serialmente, com o único custo sendo um tempo maior de processamento.

all_cores <- parallel::detectCores(logical = FALSE) - 1
registerDoFuture()

cl <- makeClusterPSOCK(all_cores)
plan(future::cluster, workers = cl)

Finalmente podemos iniciar o trabalho de tuning dos hiperparâmetros utilizando tune_grid(). A função irá estimar os 32 modelos e calcular um conjunto de medidas de performance como a acurácia e o RMSE para cada um deles. A função precisa do novo workflow (rf_wflow_tune) com o modelo e a receita a ser utilizada; dos valores possíveis dos hiperparâmetros (rf_grid) e do objeto indicando quais são os folds do cross-validation (cv_folds). Como os hiperparâmetros possuem o valor tune(), tune_grid sabe que pode atribuir valores definidos por grid_regular aos hiperparâmetros.

rf_grid_search <- tune_grid(workflow_rf_tuning, 
                            grid = rf_grid, 
                            resamples = cv_folds,
                       param_info = rf_parametros)
rf_grid_search

## # Tuning results
## # 4-fold cross-validation 
## # A tibble: 4 x 4
##   splits              id    .metrics          .notes          
##   <list>              <chr> <list>            <list>          
## 1 <split [5282/1761]> Fold1 <tibble [16 x 7]> <tibble [9 x 1]>
## 2 <split [5282/1761]> Fold2 <tibble [16 x 7]> <tibble [9 x 1]>
## 3 <split [5282/1761]> Fold3 <tibble [16 x 7]> <tibble [9 x 1]>
## 4 <split [5283/1760]> Fold4 <tibble [16 x 7]> <tibble [9 x 1]>

Alguns (bons) minutos depois, com os 32 modelos estimados, podemos calcular o valor médio de cada medida de performance para as 8 combinações diferentes de hiperparâmetros. A tabela abaixo lista os cinco modelos com melhores performance segundo o critério de ROC-UAC a partir do uso da função show_best().

rf_grid_search %>% 
  show_best(n = 5, metric = 'roc_auc') %>% 
  kable()

O modelo com menor ROC-AUC foi o Random Forest com mtry = 30, trees = 2000 e min_n = 40. Este vai ser o modelo que utilizaremos para realizar nossa previsão final.

Previsão Final

Podemos utilizar a função select_best para retornar o conjunto de hiperparâmetros do melhor modelo, como definido acima.

rf_param_final <- select_best(rf_grid_search, degree, metric = 'roc_auc')
rf_param_final %>% kable()

Usamos finalize_workflow() para criar um workflow atualizado com os valores de hiperparâmetros do modelo acima. Isso equivale a construir um novo workflow com mtry = 30, trees = 2000 e min_n = 40. Por fim, podemos ajustar o modelo utilizando toda a base de treinamento. Os folds do cross-validation já realizaram seu trabalho.

rf_wflow_final_fit <- workflow_rf_tuning %>% 
  finalize_workflow(rf_param_final) %>% 
  fit(training(tbl_treinamento_teste))
rf_wflow_final_fit

## == Workflow [trained] ==========================================================
## Preprocessor: Recipe
## Model: rand_forest()
## 
## -- Preprocessor ----------------------------------------------------------------
## 3 Recipe Steps
## 
## * step_impute_linear()
## * step_dummy()
## * step_normalize()
## 
## -- Model -----------------------------------------------------------------------
## Ranger result
## 
## Call:
##  ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~30L,      x), num.trees = ~2000L, min.node.size = min_rows(~40L, x),      num.threads = 1, verbose = FALSE, seed = sample.int(10^5,          1), probability = TRUE) 
## 
## Type:                             Probability estimation 
## Number of trees:                  2000 
## Sample size:                      5635 
## Number of independent variables:  30 
## Mtry:                             30 
## Target node size:                 40 
## Variable importance mode:         none 
## Splitrule:                        gini 
## OOB prediction error (Brier s.):  0.1366449

E realizar a previsão do modelo, calculando todas as medidas de performance.

rf_wflow_final_fit %>% 
  predict(new_data = testing(tbl_treinamento_teste)) %>% 
  bind_cols(testing(tbl_treinamento_teste) %>% select(churn)) %>% 
  mutate_all(as.factor) %>% 
  conf_mat(churn, .pred_class) %>% 
  summary() %>% 
  select(-.estimator) %>% 
  filter(.metric %in% c('precision', 'recall', 'f_meas',
                        'accuracy', 'spec', 'sens')) %>% 
  kable()

É possível observar que o modelo de random forest com parâmetros afinado não excedeu em muito o que já havia sido obtido no nosso primeiro modelo da parte 1.

Conclusão

Existe muito espaço para melhoramento do modelo acima, sobretudo em relação ao processo de feature engineering aplicado com o uso de recipes, ou mesmo com a aplicação de outros modelos alternativos.

Para além do que foi desenvolvido na parte 1 e 2 desta série, o tidymodels oferece soluções para compararação de diferentes modelos e pré-processamentos (receitas) utilizando o pacote workflowsets. Com uma lista de modelos e uma lista de diferentes receitas é possível estimar uma grande combinação de diferentes especificações, inclusive permitindo que modelos com hiperparâmetros sejam tunados.

O tidymodels também permite a criação de assembly (conjunto de modelos), seja por stacks com o uso do pacote stacks ou por bagging, com o pacote baguette, de modo a construir classificadores e previsores melhores a partir da união de diferentes algorítimos.

Na parte 3 vamos analisar o uso do workflowsetse do pacote baguette.