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 AltoO 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 rowsApó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 rowslogistic_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 rowsNote 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 rowsAssim, 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 logisticAlé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: 620ms
##
## Call:
## kknn::train.kknn(formula = quality ~ ., data = data, ks = 5)
##
## Type of response variable: nominal
## Minimal misclassification: 0.2066854
## Best kernel: optimal
## Best k: 5logistic_fit <- logistic_spec %>%
fit(quality ~ .,
data = juice(dados_recipe))
logistic_fit## parsnip model object
##
## Fit time: 31ms
##
## 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: 3998Assim, 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.806Para 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!
Tiago Mendonça
My interests include data science, machine learning, programming, photography, traveling and music.