install.packages("tidymodels")Tidymodels
O tydymodels é um meta pacote contendo 8 pacotes focados na modelagem estatística.
São eles:
- rsample: para dividr e reamostrar dados.
- parsnip: para construir modelos.
- recipes: para pré-processamento de dados.
- workflows: para agrupar o pre e pós-processamento junto do modelo.
- tune:para a otimização de hiperparâmetros e pré-processamento.
- yardstick: para medir a efetividade dos modelos.
- broom: para converter os resultados mais usuário-amigáveis.
- dials: para criar e gerir hiperparâmetros.
Como ele funciona dentro da estrutura do tidyverse, pode ser utilizado junto do pipe |> (ou %>%).
Há outros pacotes complementares que não são carregados automaticamente com o library(tidymodels).
library(tidymodels)── Attaching packages ────────────────────────────────────── tidymodels 1.1.1 ──✔ broom 1.0.5 ✔ recipes 1.0.9
✔ dials 1.2.0 ✔ rsample 1.2.0
✔ dplyr 1.1.4 ✔ tibble 3.2.1
✔ ggplot2 3.5.1 ✔ tidyr 1.3.1
✔ infer 1.0.6 ✔ tune 1.1.2
✔ modeldata 1.3.0 ✔ workflows 1.1.3
✔ parsnip 1.1.1 ✔ workflowsets 1.0.1
✔ purrr 1.0.2 ✔ yardstick 1.3.0Warning: package 'ggplot2' was built under R version 4.3.3── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ purrr::discard() masks scales::discard()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ recipes::step() masks stats::step()
• Use tidymodels_prefer() to resolve common conflicts.Modelo de regressão logística
Possui variáveis respostas categóricas. O objetivo dos modelos é prever ou a probabilidade de chegar em uma determinada categoria ou a própria categoria em si.
O conjunto de dados de exemplo é quanto a medidas da flor de três espécies de plantas do gênero Iris.
data(iris)
glimpse(iris)Rows: 150
Columns: 5
$ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9, 5.4, 4.…
$ Sepal.Width <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1, 3.7, 3.…
$ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5, 1.…
$ Petal.Width <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0.2, 0.…
$ Species <fct> setosa, setosa, setosa, setosa, setosa, setosa, setosa, s…Para modelos de classificação é interessante ver se há um equilíbrio na proporção de classes do conjunto de dados.
iris |>
count(Species) |>
mutate(prop = n/sum(n))# A tibble: 3 × 3
Species n prop
<fct> <int> <dbl>
1 setosa 50 0.333
2 versicolor 50 0.333
3 virginica 50 0.3331ª ETAPA: Corte
Separar o conjunto de dados para treinamento e teste do modelo. Como a variável Species é bem balanceada, isso pode ser feito assim:
data_split <- initial_split(iris, strata = Species)
train_data <- training(data_split)
test_data <- testing(data_split)Se contarmos a proporção de espécies nos dois novos conjuntos de dados, elas seguirão similares ao conjunto inicial.
train_data |>
count(Species) |>
mutate(prop = n/sum(n))
#> # A tibble: 3 × 3
#> Species n prop
#> <fct> <int> <dbl>
#> 1 setosa 37 0.333
#> 2 versicolor 37 0.333
#> 3 virginica 37 0.333
test_data |>
count(Species) |>
mutate(prop = n/sum(n))
#> # A tibble: 3 × 3
#> Species n prop
#> <fct> <int> <dbl>
#> 1 setosa 13 0.333
#> 2 versicolor 13 0.333
#> 3 virginica 13 0.3332ª ETAPA: Modelo de regressão
Como tem-se 3 níveis na variável resposta, um modelo de regressão logístico multinominal deve ser usado. Para criar ele será usado o pacote parsnip.
mlr_mod <-
multinom_reg(mode = "classification",
mixture = 1) |> # Para escolha de um modelo mais simples
set_engine("nnet")3ª ETAPA: Receita
Uma receita dentro do tidymodelsnada mais é do que uma lista de etapas a ser realizada para o pré-processamento dos dados antes de rodar o modelo. Isso inclui a normalização ou padronização dos dados, transformação de variáveis (ex.:, datas em meses ou dias da semana, variáveis categóricas em dummy).
iris_recipe <- recipe(Species ~ ., data = train_data) |>
step_normalize(all_numeric()) |> # Normalização
step_dummy(all_nominal(), -all_outcomes()) |> # Converção em dummy
step_zv(all_predictors()) # Remove variáveis irrelevantes4ª ETAPA: Workflow
Isso permite que o pré-processmento dos dados seja feito de uma vez dado o modelo.
logit_workflow <- workflow() |>
add_recipe(iris_recipe) |>
add_model(mlr_mod)5ª ETAPA: Treinamento
Vamos utilizar os dados de treinamento no modelo.
lr_res <- logit_workflow |>
fit(data = train_data)
lr_res|>
extract_fit_parsnip() parsnip model object
Call:
nnet::multinom(formula = ..y ~ ., data = data, trace = FALSE)
Coefficients:
(Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
versicolor 8.204878 0.3527094 -3.143007 9.689908 6.177414
virginica -8.818932 -1.5526570 -6.044161 24.144414 17.458209
Residual Deviance: 10.59006
AIC: 30.59006 6ª ETAPA: Predição e avaliação
Agora é necessário ver se o modelo se ajusta ao restante dos dados.
probabilities <- predict(lr_res, test_data) |>
bind_cols(predict(lr_res, test_data, type = "prob")) |>
bind_cols(test_data |>
select(Species))Para avaliar a performance, pode-se usar medidas de performance usando o pacote yardstick.
probabilities |>
roc_auc(truth = Species, .pred_setosa:.pred_virginica) # Melhor para avaliar a predição de probabilidade# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 roc_auc hand_till 1probabilities |>
accuracy(truth = Species, estimate = .pred_class) # Melhor para avaliar a predição de classe# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 1O modelo está muito bem ajustado para ser verdade. É necessário penalizar o modelo para que ele não cometa erros com outros conjuntos de dados semelhantes.
7ª ETAPA: Sintonização
Novos objetos
Vamos criar um outro modelo, mas agora com um argumento de penalização. Como não sabemos o quanto devemos aplicar nessa penalização, isso é considerado um hiperparâmetro a ser refinado depois.
lr_mod_p <-
multinom_reg(penalty = tune(),
mixture = 1) |>
set_mode("classification") |>
set_engine("nnet")Vamos usar a mesma receita, mas precisamos de um novo workflow:
lr_wf_p <-
workflow() |>
add_model(lr_mod_p) |>
add_recipe(iris_recipe)Gerar os valores aleatórios a serem usados no teste de penalização.
lr_reg_grid <- tibble(penalty = 10^seq(-4, -1, length.out = 30))Reamostragem
Ao invés de treinar esses modelos com o novo parâmetro no conjunto todo, vamos utilizar uma pequena parcela dos dados de treinamento.
val_set <- validation_split(train_data,
strata = Species,
prop = 0.8)Treinamento com o hiperparâmetro
Cada valor gerado anteriormente será usado no treinamento, depois é avaliado qual deles performa melhor.
lr_res_p <-
lr_wf_p |>
tune_grid(val_set,
grid = lr_reg_grid,
control = control_grid(save_pred = T), # Salva as predições do conjunto de validação
metrics = metric_set(roc_auc)) # Para avaliar a performance de cada valor de penalizaçãoWarning: package 'nnet' was built under R version 4.3.3Agora é possível plotar o resultado da predição de valores do parãmetro.
lr_plot <-
lr_res_p %>%
collect_metrics() %>%
ggplot(aes(x = penalty, y = mean)) +
geom_point() +
geom_line() +
ylab("Area under the ROC Curve") +
scale_x_log10(labels = scales::label_number())
lr_plot 
Para uma melhor visualização do valor adequado, vamos colocar isso numa tabela:
top_models <- lr_res_p |>
show_best(metric = "roc_auc", n = 15) |>
arrange(penalty)
top_models# A tibble: 15 × 7
penalty .metric .estimator mean n std_err .config
<dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
1 0.00356 roc_auc hand_till 0.990 1 NA Preprocessor1_Model16
2 0.00452 roc_auc hand_till 0.995 1 NA Preprocessor1_Model17
3 0.00574 roc_auc hand_till 0.995 1 NA Preprocessor1_Model18
4 0.00728 roc_auc hand_till 0.995 1 NA Preprocessor1_Model19
5 0.00924 roc_auc hand_till 0.995 1 NA Preprocessor1_Model20
6 0.0117 roc_auc hand_till 0.995 1 NA Preprocessor1_Model21
7 0.0149 roc_auc hand_till 0.995 1 NA Preprocessor1_Model22
8 0.0189 roc_auc hand_till 0.995 1 NA Preprocessor1_Model23
9 0.0240 roc_auc hand_till 0.995 1 NA Preprocessor1_Model24
10 0.0304 roc_auc hand_till 0.995 1 NA Preprocessor1_Model25
11 0.0386 roc_auc hand_till 0.995 1 NA Preprocessor1_Model26
12 0.0489 roc_auc hand_till 0.995 1 NA Preprocessor1_Model27
13 0.0621 roc_auc hand_till 0.990 1 NA Preprocessor1_Model28
14 0.0788 roc_auc hand_till 0.990 1 NA Preprocessor1_Model29
15 0.1 roc_auc hand_till 0.990 1 NA Preprocessor1_Model30O melhor valor foi do modelo 14 com valor de penalidade de 0.00221. No entanto, múltiplos valores de penalidade retornam um alto valor de ROC AUC.
lr_auc <-
lr_res_p |>
collect_predictions(parameters = lr_res_p|>
collect_metrics() |>
arrange(penalty) |>
slice(14)) |>
roc_curve(truth = Species, .pred_setosa:.pred_virginica) |>
mutate(model = "LR")
autoplot(lr_auc)
O modelo é bom para prever todas as espécies de Iris.
Modelo de Árvore de Decisão (randomForest)
É mais flexível que a regressão logística. Trata-se de um modelo de conjunto normalmente composto por milhares de árvores de decisão, em que cada árvore individual vê uma versão ligeiramente diferente dos dados de treinamento e aprende uma sequência de regras de divisão para prever novos dados. Cada árvore não é linear e a agregação entre árvores torna o modelo geral mais robusto e estável em comparação com as árvores individuais.
Serão usados os mesmos dados divididos anteriormente.
Modelo
rf_mod <-
rand_forest(mtry = tune(), min_n = tune(), trees = 1000) |>
set_engine("ranger", num.threads = parallel::detectCores()) |>
set_mode("classification")Workflow
rf_wf <-
workflow() |>
add_model(rf_mod) |>
add_recipe(iris_recipe)Treinamento e sintonização
Vamos ver o que precisa ser sintonizado nos hiperparâmetros definidos.
extract_parameter_set_dials(rf_mod)Collection of 2 parameters for tuning
identifier type object
mtry mtry nparam[?]
min_n min_n nparam[+]
Model parameters needing finalization:
# Randomly Selected Predictors ('mtry')
See `?dials::finalize` or `?dials::update.parameters` for more information.# `mtry` sets the number of predictor variables that each node in the decision tree “sees” and can learn about
# `min_n` sets the minimum n to split at any nodeVamos sintonizar o modelo.
set.seed(345)
rf_res <-
rf_wf|>
tune_grid(val_set,
grid = 25,
control = control_grid(save_pred = TRUE),
metrics = metric_set(roc_auc))i Creating pre-processing data to finalize unknown parameter: mtryWarning: package 'ranger' was built under R version 4.3.3Os melhores modelos:
# A tibble: 5 × 8
mtry min_n .metric .estimator mean n std_err .config
<int> <int> <chr> <chr> <dbl> <int> <dbl> <chr>
1 1 3 roc_auc hand_till 0.987 1 NA Preprocessor1_Model19
2 2 25 roc_auc hand_till 0.982 1 NA Preprocessor1_Model03
3 2 26 roc_auc hand_till 0.982 1 NA Preprocessor1_Model15
4 2 12 roc_auc hand_till 0.982 1 NA Preprocessor1_Model17
5 2 18 roc_auc hand_till 0.982 1 NA Preprocessor1_Model22
Todos os parâmetros escolhidos pelo modelo indicam um ótimo ajuste aos dados (ROC AUC > 0.9).
Predição
rf_auc <- rf_res |>
collect_predictions(parameters = rf_res |>
select_best(metric = "roc_auc")) |> # Só possível porque o modelo tinha especificado `control_grid(save_pred = TRUE)`
roc_curve(Species, .pred_setosa:.pred_virginica) |>
mutate(model = "RF")
rf_auc# A tibble: 74 × 5
.level .threshold specificity sensitivity model
<chr> <dbl> <dbl> <dbl> <chr>
1 setosa -Inf 0 1 RF
2 setosa 0 0 1 RF
3 setosa 0.000154 0.188 1 RF
4 setosa 0.0005 0.25 1 RF
5 setosa 0.000833 0.375 1 RF
6 setosa 0.00282 0.438 1 RF
7 setosa 0.00490 0.5 1 RF
8 setosa 0.005 0.562 1 RF
9 setosa 0.00649 0.625 1 RF
10 setosa 0.011 0.688 1 RF
# ℹ 64 more rowsComparando a performance dos dois modelos:
Ambos modelos performaram bem na predição das espécies.
Melhor modelo
Usando o modelo de randomForest podemos ver quais variáveis foram mais importantes na classificação.
Vamos construir um novo modelo com os parâmetros encontrados anteriormente.
last_rf_mod <-
rand_forest(mtry = 2, min_n = 7, trees = 1000) |>
set_engine("ranger", importance = "impurity", num.threads = parallel::detectCores()) |>
set_mode("classification")
last_rf_workflow <-
rf_wf |>
update_model(last_rf_mod)
set.seed(345)
last_rf_fit <-
last_rf_workflow |>
last_fit(data_split)
last_rf_fit |>
collect_metrics()# A tibble: 2 × 4
.metric .estimator .estimate .config
<chr> <chr> <dbl> <chr>
1 accuracy multiclass 0.974 Preprocessor1_Model1
2 roc_auc hand_till 1 Preprocessor1_Model1A estimativa de ROC AUC foi tão boa quanto pela regressão linear.
last_rf_fit |>
extract_fit_parsnip() |>
vip()
Assim viu-se que as medidas de pétala são mais importantes na predição das espécies que as medidas de sépalas.
