Introdução ao tidymodels

Introdução ao tidymodels

Escrevi esse post como parte de um processo para entender as possibilidades de modelagem com a abordagem utilizada no tidymodels. As principais referências para esse post são dadas pelo vídeo da Julia Silge e uma apresentação do Max Kuhn.

Vamos comparar os desempenhos preditivos dos modelos KNN e regressão logística utilizando os dados Wine Quality Data Set. Esses dados apresentam diversas medidas obtidas para os vinhos além de escores de qualidade.

Primeiro, vamos carregar os pacotes, fazer a leitura e recodificação das variáveis.

library(tidymodels)
library(tidyverse)
library(tune)

url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-white.csv"

dados <- read.csv2(url, dec = ".")

dados$quality <- ifelse(dados$quality >= 6 ,"Alto", "Baixo")

head(dados)

  fixed.acidity volatile.acidity citric.acid residual.sugar chlorides
1           7.0             0.27        0.36           20.7     0.045
2           6.3             0.30        0.34            1.6     0.049
3           8.1             0.28        0.40            6.9     0.050
4           7.2             0.23        0.32            8.5     0.058
5           7.2             0.23        0.32            8.5     0.058
6           8.1             0.28        0.40            6.9     0.050
  free.sulfur.dioxide total.sulfur.dioxide density   pH sulphates alcohol
1                  45                  170  1.0010 3.00      0.45     8.8
2                  14                  132  0.9940 3.30      0.49     9.5
3                  30                   97  0.9951 3.26      0.44    10.1
4                  47                  186  0.9956 3.19      0.40     9.9
5                  47                  186  0.9956 3.19      0.40     9.9
6                  30                   97  0.9951 3.26      0.44    10.1
  quality
1    Alto
2    Alto
3    Alto
4    Alto
5    Alto
6    Alto

O primeiro passo é separar os dados em conjunto de treinamento e de teste. Para isso, utilizaremos a função initial_split do pacote rsample. O primeiro argumento dessa funçao será dado pelo banco dados que particionaremos e o segundo argumento (prop) indicará a proporção dos dados que será utilizada para o treinamento.

set.seed(123)
dados_split <- initial_split(dados, prop = 0.8)
dados_split

<Training/Testing/Total>
<3918/980/4898>

Note que, ao imprimir o objeto dados_split, recebemos uma indicação de como foi feita a partição. Assim, temos 3.919 observações para treinamento, 979 para teste e o total de observações é dado por 4.898 vinhos. Portanto, podemos definir explicitamente os conjuntos de treinamento e teste da seguinte forma:

dados_tr   <- training(dados_split)
dados_test <- testing(dados_split)

Cooking time!

O fluxo utilizado com o pacote recipes pode ser definido pela seguinte figura (retirada da apresentação do Max Kuhn).

Primeiro definimos a receita/recipe, ou seja, como será dado o processamento dos dados, depois fazemos a preparação/prepare que especifica as estimavas necessárias para cada passo (por exemplo, para padronização ou alguma transformação como Box-Cox) e, por fim, aplicamos a receita com as estimativas utilizando as funções bake ou juice.

Agora podemos preparar a receita com o conjunto de treinamento definido. Uma receita define uma série de transformações que aplicaremos aos dados. Nesse caso criaremos uma receita de modelo com a variável quality como resposta em função das demais variáveis. Após a definição da fórmula do modelo, faremos a padronização de todas as variáveis numéricas com a função step_normalize. Além da normalização, podemos utilizar diversas tranformações. A seguir listamos alguns exemplos:

• step_BoxCox: aplica transformação de Box-Cox

• step_discretize: converte os dados númericos em fatores com níveis de tamanhos amostrais aproximadamente iguais

• step_dummy: converte variáveis nominais ou de fatores em variáveis indicadoras

• step_log: aplica logaritmo

• step_meaninput, step_medianinput, step_modeinput: imputa os dados faltantes por uma das medidas (média, mediana e moda) do conjunto de treinamento

• step_normalize: padroniza as variáveis para ter média zero e desvio padrão igual a um

• step_YeoJohnson: transforma as variáveis de acordo com a transformação de Yeo-Johnson

Após definir a fórmula utilizada e a padronização dos dados, utilizaremos a função prep para aplicar as transformações a um determinado conjunto de dados. Nesse caso, aplicaremos no conjunto de treinamento (dados_tr).

dados_recipe <- dados_tr %>% 
  recipe(quality ~ .) %>%
  step_normalize(all_numeric()) %>%
  prep()

dados_recipe

Recipe

Inputs:

      role #variables
   outcome          1
 predictor         11

Training data contained 3918 data points and no missing data.

Operations:

Centering and scaling for fixed.acidity, volatile.acidity, citric.acid, r... [trained]

Para obter os dados após a aplicação da receita, podemos utilizar os seguintes comandos:

juice(dados_recipe)

# A tibble: 3,918 × 12
   fixed.acid…¹ volat…² citri…³ resid…⁴ chlor…⁵ free.…⁶ total…⁷ density       pH
          <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>    <dbl>
 1        1.74  -0.285  -0.0322 -0.756   0.330  -1.35   -1.57    0.0166 -1.97   
 2       -0.770  0.0111  0.956   0.235  -0.0335  0.620   1.53    0.581   0.463  
 3        1.02   0.0111  2.03    1.15    0.0119  1.43    0.915   1.44   -0.719  
 4       -1.13   0.405   0.462  -0.598   3.10   -0.653   0.0404 -0.197  -0.654  
 5       -1.37  -0.186  -0.773   2.28    0.602  -0.711  -0.338   1.25    0.726  
 6       -0.173 -0.975  -0.444   0.770  -0.306  -0.364  -0.550   0.234  -0.522  
 7        0.425 -0.383  -0.362  -0.816  -0.397   0.0990 -0.857  -0.665   0.529  
 8        2.93   0.603   1.29    0.195   0.284  -0.306   0.182   1.34   -0.128  
 9       -1.01   0.110  -0.691  -0.974  -0.578  -0.306  -0.574  -2.02   -0.719  
10       -0.292 -1.27   -0.115  -0.0625 -0.578   1.37   -0.243  -0.692   0.00315
# … with 3,908 more rows, 3 more variables: sulphates <dbl>, alcohol <dbl>,
#   quality <fct>, and abbreviated variable names ¹​fixed.acidity,
#   ²​volatile.acidity, ³​citric.acid, ⁴​residual.sugar, ⁵​chlorides,
#   ⁶​free.sulfur.dioxide, ⁷​total.sulfur.dioxide

bake(dados_recipe, new_data = dados_tr)

# A tibble: 3,918 × 12
   fixed.acid…¹ volat…² citri…³ resid…⁴ chlor…⁵ free.…⁶ total…⁷ density       pH
          <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>    <dbl>
 1        1.74  -0.285  -0.0322 -0.756   0.330  -1.35   -1.57    0.0166 -1.97   
 2       -0.770  0.0111  0.956   0.235  -0.0335  0.620   1.53    0.581   0.463  
 3        1.02   0.0111  2.03    1.15    0.0119  1.43    0.915   1.44   -0.719  
 4       -1.13   0.405   0.462  -0.598   3.10   -0.653   0.0404 -0.197  -0.654  
 5       -1.37  -0.186  -0.773   2.28    0.602  -0.711  -0.338   1.25    0.726  
 6       -0.173 -0.975  -0.444   0.770  -0.306  -0.364  -0.550   0.234  -0.522  
 7        0.425 -0.383  -0.362  -0.816  -0.397   0.0990 -0.857  -0.665   0.529  
 8        2.93   0.603   1.29    0.195   0.284  -0.306   0.182   1.34   -0.128  
 9       -1.01   0.110  -0.691  -0.974  -0.578  -0.306  -0.574  -2.02   -0.719  
10       -0.292 -1.27   -0.115  -0.0625 -0.578   1.37   -0.243  -0.692   0.00315
# … with 3,908 more rows, 3 more variables: sulphates <dbl>, alcohol <dbl>,
#   quality <fct>, and abbreviated variable names ¹​fixed.acidity,
#   ²​volatile.acidity, ³​citric.acid, ⁴​residual.sugar, ⁵​chlorides,
#   ⁶​free.sulfur.dioxide, ⁷​total.sulfur.dioxide

Agora podemos aplicar essa receita aos dados de teste que será utilizado posteriormente para avaliação do desempenho preditivo.

test_proc <- dados_recipe %>% 
              bake(new_data = dados_test)

A seguir vamos criar a especificação dos modelos preditivos.

knn_spec <- nearest_neighbor() %>% 
             set_engine("kknn") %>% 
             set_mode("classification")

logistic_spec <- logistic_reg() %>% 
                  set_engine("glm") %>% 
                  set_mode("classification")

A função set_engine especifica qual pacote ou sistema será utilizado para ajustar o modelo. Aqui utilizaremos o kknn e glm. Além disso, o modo do modelo é classificação. Poderíamos utilizar o modo regression, caso a variável resposta fosse numérica.

Uma vez estabelecida a especificação do modelo, podemos utilizar a função mc_cv (Monte Carlo Cross-Validation) para fazer várias amostras de treinamento e teste com base no conjunto de dados de treinamento processado. O argumento prop indica a proporção utilizada para modelagem. O padrão da função é repetir esse processo 25 vezes. O resultado é um tibble com as indicações do número de observações para treinamento, teste e total na coluna splits e uma coluna id atribuindo um código para cada reamostragem.

set.seed(1234)

validation_splits <- mc_cv(juice(dados_recipe), prop = 0.8)
validation_splits

# Monte Carlo cross-validation (0.8/0.2) with 25 resamples  
# A tibble: 25 × 2
   splits             id        
   <list>             <chr>     
 1 <split [3134/784]> Resample01
 2 <split [3134/784]> Resample02
 3 <split [3134/784]> Resample03
 4 <split [3134/784]> Resample04
 5 <split [3134/784]> Resample05
 6 <split [3134/784]> Resample06
 7 <split [3134/784]> Resample07
 8 <split [3134/784]> Resample08
 9 <split [3134/784]> Resample09
10 <split [3134/784]> Resample10
# … with 15 more rows

Após a criação de várias amostras de treinamento e teste no objeto validation_splits, podemos utilizar a função fit_resamples do pacote tune. Com isso, é possível ajustar os modelos knn e logístico com as especificações guardadas nos objetos knn_spec e logistic_spec de acordo com a seguinte função:

knn_res <- fit_resamples(knn_spec, 
                         quality ~ ., 
                         validation_splits, 
                         control = control_resamples(save_pred = TRUE))

knn_res

# Resampling results
# Monte Carlo cross-validation (0.8/0.2) with 25 resamples  
# A tibble: 25 × 5
   splits             id         .metrics         .notes           .predictions
   <list>             <chr>      <list>           <list>           <list>      
 1 <split [3134/784]> Resample01 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 2 <split [3134/784]> Resample02 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 3 <split [3134/784]> Resample03 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 4 <split [3134/784]> Resample04 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 5 <split [3134/784]> Resample05 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 6 <split [3134/784]> Resample06 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 7 <split [3134/784]> Resample07 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 8 <split [3134/784]> Resample08 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 9 <split [3134/784]> Resample09 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
10 <split [3134/784]> Resample10 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
# … with 15 more rows

logistic_res <- fit_resamples(logistic_spec,  
                              quality ~ ., 
                              validation_splits, 
                              control = control_resamples(save_pred = TRUE))

logistic_res

# Resampling results
# Monte Carlo cross-validation (0.8/0.2) with 25 resamples  
# A tibble: 25 × 5
   splits             id         .metrics         .notes           .predictions
   <list>             <chr>      <list>           <list>           <list>      
 1 <split [3134/784]> Resample01 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 2 <split [3134/784]> Resample02 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 3 <split [3134/784]> Resample03 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 4 <split [3134/784]> Resample04 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 5 <split [3134/784]> Resample05 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 6 <split [3134/784]> Resample06 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 7 <split [3134/784]> Resample07 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 8 <split [3134/784]> Resample08 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
 9 <split [3134/784]> Resample09 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
10 <split [3134/784]> Resample10 <tibble [2 × 4]> <tibble [0 × 3]> <tibble>    
# … with 15 more rows

Note que agora, além das colunas splits e id, temos as colunas .metrics, .notes e .predictions. A coluna .metrics apresenta as medidas de desempenho calculadas para cada conjunto gerado e a coluna .predictions guarda a probabilidade predita para cada classe e a classificação observada (as colunas com as probabilidades de cada classe serão nomeadas de acordo com o padrão .pred_classe)

knn_res %>% 
  select(id, .metrics) %>% 
  unnest(.metrics)

# A tibble: 50 × 5
   id         .metric  .estimator .estimate .config             
   <chr>      <chr>    <chr>          <dbl> <chr>               
 1 Resample01 accuracy binary         0.800 Preprocessor1_Model1
 2 Resample01 roc_auc  binary         0.850 Preprocessor1_Model1
 3 Resample02 accuracy binary         0.793 Preprocessor1_Model1
 4 Resample02 roc_auc  binary         0.819 Preprocessor1_Model1
 5 Resample03 accuracy binary         0.811 Preprocessor1_Model1
 6 Resample03 roc_auc  binary         0.845 Preprocessor1_Model1
 7 Resample04 accuracy binary         0.797 Preprocessor1_Model1
 8 Resample04 roc_auc  binary         0.848 Preprocessor1_Model1
 9 Resample05 accuracy binary         0.788 Preprocessor1_Model1
10 Resample05 roc_auc  binary         0.830 Preprocessor1_Model1
# … with 40 more rows

Assim, pode ser interessante verificar visualmente o desempenho preditivo de cada modelo gerado. Para isso, podemos utilizar o seguinte código:

knn_res %>% 
  select(id, .metrics) %>% 
  unnest(.metrics) %>% 
  mutate(model = "knn") %>% 
  bind_rows(logistic_res %>% 
            select(id, .metrics) %>% 
            unnest(.metrics) %>% 
            mutate(model = "logistic")) %>% 
  ggplot(aes(id, .estimate, group = model, color = model)) + 
    geom_point(size = 1.2) + 
    facet_wrap(~.metric) + 
    coord_flip()

Podemos obter, de forma mais prática, um resumo dessas observações com a função collect_metrics.

knn_res %>% 
  collect_metrics() %>% 
  mutate(model = "knn") %>% 
  bind_rows(logistic_res %>% 
            collect_metrics() %>% 
            mutate(model = "logistic"))

# A tibble: 4 × 7
  .metric  .estimator  mean     n std_err .config              model   
  <chr>    <chr>      <dbl> <int>   <dbl> <chr>                <chr>   
1 accuracy binary     0.789    25 0.00304 Preprocessor1_Model1 knn     
2 roc_auc  binary     0.827    25 0.00327 Preprocessor1_Model1 knn     
3 accuracy binary     0.748    25 0.00336 Preprocessor1_Model1 logistic
4 roc_auc  binary     0.794    25 0.00341 Preprocessor1_Model1 logistic

Além da medida resumo, podemos utilizar as previsões obtidas para os 25 conjuntos de teste gerados para fazer uma curva ROC. Para isso, usamos o comando unnest, definimos a curva roc com o comando roc_curve considerando o valor observado como primeiro argumento e a probabilidade do vinho apresentar alta qualidade no segundo argumento. Após esses passos, utilizamos a função autoplot.

knn_res %>% 
  unnest(.predictions) %>% 
  mutate(model = "knn") %>% 
  bind_rows(logistic_res %>% 
            unnest(.predictions) %>% 
            mutate(model = "logistic")) %>% 
  group_by(model) %>% 
  roc_curve(quality, .pred_Alto) %>% 
  autoplot()

Após o estudo dos modelos com os conjuntos de treinamentos obtidos por reamostragem, podemos ajustá-los, agora considerando o conjunto de treinamento completo, com a especificação dada por knn_spec/logistic_spec. Assim, teremos

knn_fit <- knn_spec %>% 
            fit(quality ~ ., 
                data = juice(dados_recipe))

knn_fit

parsnip model object


Call:
kknn::train.kknn(formula = quality ~ ., data = data, ks = min_rows(5,     data, 5))

Type of response variable: nominal
Minimal misclassification: 0.1993364
Best kernel: optimal
Best k: 5

logistic_fit <- logistic_spec %>% 
                  fit(quality ~ ., 
                  data = juice(dados_recipe))

logistic_fit

parsnip model object


Call:  stats::glm(formula = quality ~ ., family = stats::binomial, data = data)

Coefficients:
         (Intercept)         fixed.acidity      volatile.acidity  
           -0.915911              0.013944              0.647401  
         citric.acid        residual.sugar             chlorides  
           -0.009177             -0.730209             -0.014389  
 free.sulfur.dioxide  total.sulfur.dioxide               density  
           -0.098973              0.013988              0.607163  
                  pH             sulphates               alcohol  
           -0.113490             -0.180647             -1.030140  

Degrees of Freedom: 3917 Total (i.e. Null);  3906 Residual
Null Deviance:      5000 
Residual Deviance: 3964     AIC: 3988

Assim, para avaliar a AUC obtida com o conjunto de teste, utilizamos:

knn_fit %>% 
  predict(new_data = test_proc, type = "prob") %>% 
  mutate(truth = test_proc$quality, 
         model = "knn") %>%
  bind_rows(logistic_fit %>% 
            predict(new_data = test_proc, type = "prob") %>% 
            mutate(truth = test_proc$quality, 
                   model = "logistic")) %>% 
  group_by(model) %>% 
  roc_auc(truth, .pred_Alto) 

# A tibble: 2 × 4
  model    .metric .estimator .estimate
  <chr>    <chr>   <chr>          <dbl>
1 knn      roc_auc binary         0.833
2 logistic roc_auc binary         0.805

Para fazer a curva ROC, poderíamos utilizar diretamente:

knn_fit %>% 
  predict(new_data = test_proc, type = "prob") %>% 
  mutate(truth = test_proc$quality, 
         model = "knn") %>%
  bind_rows(logistic_fit %>% 
            predict(new_data = test_proc, type = "prob") %>% 
            mutate(truth = test_proc$quality, 
                   model = "logistic")) %>% 
  group_by(model) %>% 
  roc_curve(truth, .pred_Alto) %>% 
  autoplot()

O fluxo apresentado acima pode não ser exatamente o fluxo utilizado num processo de modelagem, mas apresenta as principais formas de aplicação dessas abordagens no contexto do tidymodels.

Caso tenha alguma crítica, sugestão ou comentário, me envie uma mensagem!