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
## <Analysis/Assess/Total>
## <3919/979/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-Coxstep_discretize
: converte os dados númericos em fatores com níveis de tamanhos amostrais aproximadamente iguaisstep_dummy
: converte variáveis nominais ou de fatores em variáveis indicadorasstep_log
: aplica logaritmostep_meaninput
,step_medianinput
,step_modeinput
: imputa os dados faltantes por uma das medidas (média, mediana e moda) do conjunto de treinamentostep_normalize
: padroniza as variáveis para ter média zero e desvio padrão igual a umstep_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
## Data Recipe
##
## Inputs:
##
## role #variables
## outcome 1
## predictor 11
##
## Training data contained 3919 data points and no missing data.
##
## Operations:
##
## Centering and scaling for fixed.acidity, volatile.acidity, ... [trained]
Para obter os dados após a aplicação da receita, podemos utilizar os seguintes comandos:
juice(dados_recipe)
## # A tibble: 3,919 x 12
## fixed.acidity volatile.acidity citric.acid residual.sugar chlorides
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -0.677 0.200 0.0394 -0.955 0.145
## 2 1.46 0.00482 0.530 0.0772 0.190
## 3 0.390 -0.482 -0.124 0.389 0.552
## 4 1.46 0.00482 0.530 0.0772 0.190
## 5 -0.796 0.394 -1.43 0.0966 -0.0357
## 6 0.153 -0.0925 0.203 2.76 -0.0357
## 7 -0.677 0.200 0.0394 -0.955 0.145
## 8 1.46 -0.579 0.775 -0.974 -0.0810
## 9 1.46 -0.0925 0.612 -0.984 -0.578
## 10 2.05 -0.482 0.530 -0.448 -0.488
## # ... with 3,909 more rows, and 7 more variables: free.sulfur.dioxide <dbl>,
## # total.sulfur.dioxide <dbl>, density <dbl>, pH <dbl>, sulphates <dbl>,
## # alcohol <dbl>, quality <fct>
bake(dados_recipe, new_data = dados_tr)
## # A tibble: 3,919 x 12
## fixed.acidity volatile.acidity citric.acid residual.sugar chlorides
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -0.677 0.200 0.0394 -0.955 0.145
## 2 1.46 0.00482 0.530 0.0772 0.190
## 3 0.390 -0.482 -0.124 0.389 0.552
## 4 1.46 0.00482 0.530 0.0772 0.190
## 5 -0.796 0.394 -1.43 0.0966 -0.0357
## 6 0.153 -0.0925 0.203 2.76 -0.0357
## 7 -0.677 0.200 0.0394 -0.955 0.145
## 8 1.46 -0.579 0.775 -0.974 -0.0810
## 9 1.46 -0.0925 0.612 -0.984 -0.578
## 10 2.05 -0.482 0.530 -0.448 -0.488
## # ... with 3,909 more rows, and 7 more variables: free.sulfur.dioxide <dbl>,
## # total.sulfur.dioxide <dbl>, density <dbl>, pH <dbl>, sulphates <dbl>,
## # alcohol <dbl>, quality <fct>
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 x 2
## splits id
## <list> <chr>
## 1 <split [3.1K/783]> Resample01
## 2 <split [3.1K/783]> Resample02
## 3 <split [3.1K/783]> Resample03
## 4 <split [3.1K/783]> Resample04
## 5 <split [3.1K/783]> Resample05
## 6 <split [3.1K/783]> Resample06
## 7 <split [3.1K/783]> Resample07
## 8 <split [3.1K/783]> Resample08
## 9 <split [3.1K/783]> Resample09
## 10 <split [3.1K/783]> 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(quality ~ .,
knn_spec,
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 x 5
## splits id .metrics .notes .predictions
## <list> <chr> <list> <list> <list>
## 1 <split [3.1K/783~ Resample01 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 2 <split [3.1K/783~ Resample02 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 3 <split [3.1K/783~ Resample03 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 4 <split [3.1K/783~ Resample04 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 5 <split [3.1K/783~ Resample05 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 6 <split [3.1K/783~ Resample06 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 7 <split [3.1K/783~ Resample07 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 8 <split [3.1K/783~ Resample08 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 9 <split [3.1K/783~ Resample09 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 10 <split [3.1K/783~ Resample10 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## # ... with 15 more rows
logistic_res <- fit_resamples(quality ~ .,
logistic_spec,
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 x 5
## splits id .metrics .notes .predictions
## <list> <chr> <list> <list> <list>
## 1 <split [3.1K/783~ Resample01 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 2 <split [3.1K/783~ Resample02 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 3 <split [3.1K/783~ Resample03 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 4 <split [3.1K/783~ Resample04 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 5 <split [3.1K/783~ Resample05 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 6 <split [3.1K/783~ Resample06 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 7 <split [3.1K/783~ Resample07 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 8 <split [3.1K/783~ Resample08 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 9 <split [3.1K/783~ Resample09 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## 10 <split [3.1K/783~ Resample10 <tibble [2 x 3~ <tibble [0 x ~ <tibble [783 x 5~
## # ... 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 x 4
## id .metric .estimator .estimate
## <chr> <chr> <chr> <dbl>
## 1 Resample01 accuracy binary 0.766
## 2 Resample01 roc_auc binary 0.815
## 3 Resample02 accuracy binary 0.785
## 4 Resample02 roc_auc binary 0.816
## 5 Resample03 accuracy binary 0.760
## 6 Resample03 roc_auc binary 0.804
## 7 Resample04 accuracy binary 0.778
## 8 Resample04 roc_auc binary 0.804
## 9 Resample05 accuracy binary 0.806
## 10 Resample05 roc_auc binary 0.826
## # ... 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 x 6
## .metric .estimator mean n std_err model
## <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 accuracy binary 0.783 25 0.00244 knn
## 2 roc_auc binary 0.820 25 0.00261 knn
## 3 accuracy binary 0.750 25 0.00283 logistic
## 4 roc_auc binary 0.800 25 0.00336 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
##
## Fit time: 571ms
##
## Call:
## kknn::train.kknn(formula = quality ~ ., data = data, ks = 5)
##
## Type of response variable: nominal
## Minimal misclassification: 0.2066854
## Best kernel: optimal
## Best k: 5
logistic_fit <- logistic_spec %>%
fit(quality ~ .,
data = juice(dados_recipe))
logistic_fit
## parsnip model object
##
## Fit time: 10ms
##
## Call: stats::glm(formula = quality ~ ., family = stats::binomial, data = data)
##
## Coefficients:
## (Intercept) fixed.acidity volatile.acidity
## -0.887950 0.032699 0.662796
## citric.acid residual.sugar chlorides
## 0.010212 -0.708930 0.006118
## free.sulfur.dioxide total.sulfur.dioxide density
## -0.201736 0.066582 0.559808
## pH sulphates alcohol
## -0.122876 -0.193902 -1.020903
##
## Degrees of Freedom: 3918 Total (i.e. Null); 3907 Residual
## Null Deviance: 5021
## Residual Deviance: 3974 AIC: 3998
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 x 4
## model .metric .estimator .estimate
## <chr> <chr> <chr> <dbl>
## 1 knn roc_auc binary 0.850
## 2 logistic roc_auc binary 0.806
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!