← Voltar para a página inicial

1 Introdução

Este tutorial é parte no material do minicurso ministrado durante o I Congresso do IMECC.

O código fonte deste material está disponível em https://github.com/samarafk/CongressoIMECC2025.

2 Conjunto de dados

Vamos carregar o pacote tidyverse e ler o conjunto de dados abalone que usaremos neste tutorial:

library(tidyverse)

abalone <- read_csv("abalone.csv")

Algumas informações rápidas dos dados:

abalone %>% glimpse()
## Rows: 4,177
## Columns: 9
## $ sex            <chr> "male", "male", "female", "male", "infant", "infant", "…
## $ length         <dbl> 0.455, 0.350, 0.530, 0.440, 0.330, 0.425, 0.530, 0.545,…
## $ diameter       <dbl> 0.365, 0.265, 0.420, 0.365, 0.255, 0.300, 0.415, 0.425,…
## $ height         <dbl> 0.095, 0.090, 0.135, 0.125, 0.080, 0.095, 0.150, 0.125,…
## $ whole_weight   <dbl> 0.5140, 0.2255, 0.6770, 0.5160, 0.2050, 0.3515, 0.7775,…
## $ shucked_weight <dbl> 0.2245, 0.0995, 0.2565, 0.2155, 0.0895, 0.1410, 0.2370,…
## $ viscera_weight <dbl> 0.1010, 0.0485, 0.1415, 0.1140, 0.0395, 0.0775, 0.1415,…
## $ shell_weight   <dbl> 0.150, 0.070, 0.210, 0.155, 0.055, 0.120, 0.330, 0.260,…
## $ rings          <dbl> 15, 7, 9, 10, 7, 8, 20, 16, 9, 19, 14, 10, 11, 10, 10, …

Temos \(n=4177\) observações. A variável resposta é rings.

As primeiras 10 observações:

abalone %>% head(10)
## # A tibble: 10 × 9
##    sex    length diameter height whole_weight shucked_weight viscera_weight
##    <chr>   <dbl>    <dbl>  <dbl>        <dbl>          <dbl>          <dbl>
##  1 male    0.455    0.365  0.095        0.514         0.224          0.101 
##  2 male    0.35     0.265  0.09         0.226         0.0995         0.0485
##  3 female  0.53     0.42   0.135        0.677         0.256          0.142 
##  4 male    0.44     0.365  0.125        0.516         0.216          0.114 
##  5 infant  0.33     0.255  0.08         0.205         0.0895         0.0395
##  6 infant  0.425    0.3    0.095        0.352         0.141          0.0775
##  7 female  0.53     0.415  0.15         0.778         0.237          0.142 
##  8 female  0.545    0.425  0.125        0.768         0.294          0.150 
##  9 male    0.475    0.37   0.125        0.509         0.216          0.112 
## 10 female  0.55     0.44   0.15         0.894         0.314          0.151 
## # ℹ 2 more variables: shell_weight <dbl>, rings <dbl>

3 Data Splitting

Agora vamos carregar o pacote tidymodels e separar os dados em duas partes.

Para dividir os dados em conjuntos de treinamento e teste usamos a função initial_split(). Por padrão, a função separa 75% dos dados para o treinamento e 25% para o teste. Aqui, como exemplo, estamos usando 80% para o treinamento e isso é especificado com prop = 0.8.

A amostragem estratificada é feita através do argumento strata. A estratificação garante que os dados de teste mantenham uma distribuição próxima a dos dados de treinamento.

Como a separação em treinamento e teste é feita de forma aleatorizada é importante usar set.seed() para garantirmos sempre a mesma divisão dos dados ao executarmos o código novamente:

library(tidymodels)
set.seed(1234)
ring_split <- initial_split(abalone, prop = 0.8, strata = rings)

As funções training() e testing() são usadas para extrair os conjuntos de treinamento e teste, respectivamente:

ring_train <- training(ring_split)
ring_test <- testing(ring_split)

Mais detalhes sobre divisão do conjunto de dados em treinamento e teste são discutidos nos livros de Hastie, Tibshirani, and Friedman (2001) e James et al. (2013).

4 Modelo

No tidymodels, temos alguns passos para especificar um modelo:

  1. Escolha um modelo (model)
  2. Especifique um mecanismo (engine)
  3. Defina o modo (mode)

Por exemplo, se quisermos especificar um modelo de regressão linear:

linear_reg()
## Linear Regression Model Specification (regression)
## 
## Computational engine: lm

Depois que a forma funcional do modelo é especificada, é necessário pensar em um mecanismo/método para ajustar ou treinar o modelo, chamado de engine.

args(linear_reg)
## function (mode = "regression", engine = "lm", penalty = NULL, 
##     mixture = NULL) 
## NULL

Veja que o método default para linear_reg() é lm, que ajusta o modelo por mínimos quadrados ordinários, mas é possível escolher outros métodos.

Por exemplo, um modelo de regressão linear via mínimos quadrados generalizados:

linear_reg() %>% 
   set_engine("gls") 
## Linear Regression Model Specification (regression)
## 
## Computational engine: gls

Todos os modelos disponíveis são listados no site: https://www.tidymodels.org/find/parsnip/

5 Receitas

Antes de proceder para o ajuste/treinamento do modelo, faz-se o pré-processamento dos dados, já que:

  • Alguns modelos exigem que os preditores tenham certas características ou certos formatos.

  • Algumas variáveis requerem algum tipo de transformação.

Para isso, o tidymodels usa o que é chamado de receita (recipe).

Uma primeira receita:

ring_rec <- recipe(rings ~ ., data = ring_train)

No exemplo acima, a função recipe() apenas define a variável resposta e os preditores através da fórmula.

ring_rec %>% summary()
## # A tibble: 9 × 4
##   variable       type      role      source  
##   <chr>          <list>    <chr>     <chr>   
## 1 sex            <chr [3]> predictor original
## 2 length         <chr [2]> predictor original
## 3 diameter       <chr [2]> predictor original
## 4 height         <chr [2]> predictor original
## 5 whole_weight   <chr [2]> predictor original
## 6 shucked_weight <chr [2]> predictor original
## 7 viscera_weight <chr [2]> predictor original
## 8 shell_weight   <chr [2]> predictor original
## 9 rings          <chr [2]> outcome   original

Os passos de pré-processamento de uma receita usam os dados de treinamento para os cálculos. Tipos de cálculos no processamento: * Níveis de uma variável categórica * Avaliar se uma coluna tem variância nula (constante) - step_zv() * Calcular média e desvio-padrão para a normalização * Projetar os novos dados nas componentes principais obtidas pelos dados de treinamento.

Um outro exemplo de receita:

ring_rec <- recipe(rings ~ ., data = ring_train) %>%
  step_normalize(all_numeric_predictors()) %>%
  step_corr(all_numeric_predictors(), threshold = 0.9) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_poly(shucked_weight, degree = 2)
A tabela a seguir apresenta os dados pré-processamos segundo a receita acima:
rings height shell_weight sex_infant sex_male shucked_weight_poly_1 shucked_weight_poly_2
8 -1.044 -0.855 1 0 -0.017 0.008
7 -1.278 -0.890 1 0 -0.021 0.015
8 -0.810 -0.819 0 1 -0.010 -0.003
5 -2.214 -1.573 1 0 -0.026 0.027
5 -1.980 -1.630 1 0 -0.026 0.028
4 -2.097 -1.608 1 0 -0.027 0.030
7 -1.044 -1.178 1 0 -0.021 0.017
6 -1.629 -1.393 1 0 -0.022 0.019
7 -0.927 -0.926 0 1 -0.015 0.004
8 0.009 -0.136 0 0 -0.008 -0.006

6 Workflow

Gerenciar um modelo especificado com o parsnip (model + engine + mode) e os passos de pré-processamento usando recipes, pode ser um desafio.

Para isso, o tidymodels propõe o uso da função workflow(), que possui dois argumentos opcionais:

  • preprocessor: pode ser uma fórmula ou uma receita
  • spec: modelo especificado com parsnip

Vamos especificar um modelo de regressão linear:

reg_spec <- linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression")

A receita a ser aplicada está no objeto ring_rec.

Vamos agregar essas duas informações no workflow:

reg_wf <- workflow(ring_rec, reg_spec)

A função fit() pode ser usada para ajustar o modelo usando os dados de treinamento:

reg_fit <- reg_wf %>% fit(data = ring_train)
reg_fit
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 4 Recipe Steps
## 
## • step_normalize()
## • step_corr()
## • step_dummy()
## • step_poly()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:
## stats::lm(formula = ..y ~ ., data = data)
## 
## Coefficients:
##           (Intercept)                 height           shell_weight  
##              10.27917                0.59969                3.20871  
##            sex_infant               sex_male  shucked_weight_poly_1  
##              -1.02274               -0.02678             -127.72989  
## shucked_weight_poly_2  
##              -7.59663

A função tidy() do pacote broom resume as informações mais importantes de um modelo usando o conceito tidy:

reg_fit %>% tidy(conf.int = TRUE)
## # A tibble: 7 × 7
##   term                 estimate std.error statistic   p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)           1.03e+1    0.0746   138.    0           10.1      10.4  
## 2 height                6.00e-1    0.0704     8.52  2.35e- 17    0.462     0.738
## 3 shell_weight          3.21e+0    0.0953    33.7   2.47e-214    3.02      3.40 
## 4 sex_infant           -1.02e+0    0.118     -8.65  7.83e- 18   -1.25     -0.791
## 5 sex_male             -2.68e-2    0.0964    -0.278 7.81e-  1   -0.216     0.162
## 6 shucked_weight_poly… -1.28e+2    5.29     -24.2   6.75e-119 -138.     -117.   
## 7 shucked_weight_poly… -7.60e+0    2.50      -3.04  2.40e-  3  -12.5      -2.69

A função predict() calcula os valores preditos para o conjunto especificado:

reg_fit %>% predict(ring_train)
## # A tibble: 3,340 × 1
##    .pred
##    <dbl>
##  1  8.00
##  2  8.15
##  3  8.48
##  4  5.94
##  5  5.94
##  6  6.01
##  7  7.43
##  8  6.49
##  9  8.55
## 10 10.9 
## # ℹ 3,330 more rows

A predição com tidymodels garante que:

  • as predições estarão sempre dentro de um dataframe/tibble;
  • os nomes e tipos das colunas são previsíveis e intuitivas;
  • o número de linhas em new_data e a saída são iguais.

Pode-se usar também a função augment(), que calcula os valores preditos e resíduos, adicionando-os em uma coluna no conjunto de dados especificado:

reg_fit %>% augment(ring_train)
## # A tibble: 3,340 × 11
##    .pred    .resid sex    length diameter height whole_weight shucked_weight
##    <dbl>     <dbl> <chr>   <dbl>    <dbl>  <dbl>        <dbl>          <dbl>
##  1  8.00 -0.000960 infant  0.425    0.3    0.095        0.352         0.141 
##  2  8.15 -1.15     infant  0.355    0.28   0.085        0.290         0.095 
##  3  8.48 -0.482    male    0.465    0.355  0.105        0.480         0.227 
##  4  5.94 -0.936    infant  0.24     0.175  0.045        0.07          0.0315
##  5  5.94 -0.942    infant  0.205    0.15   0.055        0.042         0.0255
##  6  6.01 -2.01     infant  0.21     0.15   0.05         0.042         0.0175
##  7  7.43 -0.429    infant  0.39     0.295  0.095        0.203         0.0875
##  8  6.49 -0.490    infant  0.325    0.245  0.07         0.161         0.0755
##  9  8.55 -1.55     male    0.405    0.31   0.1          0.385         0.173 
## 10 10.9  -2.92     female  0.5      0.4    0.14         0.662         0.256 
## # ℹ 3,330 more rows
## # ℹ 3 more variables: viscera_weight <dbl>, shell_weight <dbl>, rings <dbl>

Para visualizar as estimativas dos parâmetros do modelo ajustado:

reg_fit %>% tidy()
## # A tibble: 7 × 5
##   term                   estimate std.error statistic   p.value
##   <chr>                     <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)             10.3       0.0746   138.    0        
## 2 height                   0.600     0.0704     8.52  2.35e- 17
## 3 shell_weight             3.21      0.0953    33.7   2.47e-214
## 4 sex_infant              -1.02      0.118     -8.65  7.83e- 18
## 5 sex_male                -0.0268    0.0964    -0.278 7.81e-  1
## 6 shucked_weight_poly_1 -128.        5.29     -24.2   6.75e-119
## 7 shucked_weight_poly_2   -7.60      2.50      -3.04  2.40e-  3

7 Avaliar e comparar modelos

Até aqui, fizemos o pré-processamento, definimos e treinamos o modelo escolhido usando os dados de treinamento.

7.1 Métricas

Como avaliar se o modelo tem bom desempenho (foco em predição)?

Olhar os resultados para cada observação não é produtivo:

reg_fit %>%
  augment(ring_train) %>% 
  head()
## # A tibble: 6 × 11
##   .pred    .resid sex    length diameter height whole_weight shucked_weight
##   <dbl>     <dbl> <chr>   <dbl>    <dbl>  <dbl>        <dbl>          <dbl>
## 1  8.00 -0.000960 infant  0.425    0.3    0.095        0.352         0.141 
## 2  8.15 -1.15     infant  0.355    0.28   0.085        0.290         0.095 
## 3  8.48 -0.482    male    0.465    0.355  0.105        0.480         0.227 
## 4  5.94 -0.936    infant  0.24     0.175  0.045        0.07          0.0315
## 5  5.94 -0.942    infant  0.205    0.15   0.055        0.042         0.0255
## 6  6.01 -2.01     infant  0.21     0.15   0.05         0.042         0.0175
## # ℹ 3 more variables: viscera_weight <dbl>, shell_weight <dbl>, rings <dbl>

Temos algumas métricas para comparar os valores preditos com os observados (erro de predição):

  • Erro Quadrático Médio: \(RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^n(y_i-\hat{y}_i)^2}\)
  • Coeficiente de determinação: \(R^2\)
  • Erro Médio Absoluto: \(MAE=\frac{1}{n}\sum_{i=1}^n|y_i-\hat{y}_i|\)

Dentro de tidymodels temos a função metrics() do pacote yardstick para avaliar o modelo. Temos que especificar os seguintes argumentos em metrics():

  • truth: nome da variável com os valores observados da resposta
  • estimate: nome da variável contendo os valores preditos
reg_fit %>%
  augment(new_data = ring_train) %>%
  metrics(rings, .pred)
## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard       2.26 
## 2 rsq     standard       0.510
## 3 mae     standard       1.62

Podemos especificar apenas uma métrica:

reg_fit %>%
  augment(new_data = ring_train) %>%
  rmse(rings, .pred)
## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        2.26

Também é possível especificar um conjunto de métricas. No exemplo a seguir, abalone_metrics é definido como um conjunto de quatro métricas: RMSE, MAE, MAPE e \(R^2\):

abalone_metrics <- metric_set(rmse, mae, mape, rsq)

E podemos avaliar este conjunto de métricas no modelo ajustado:

augment(reg_fit, new_data = ring_train) %>%
  abalone_metrics(rings, .pred)
## # A tibble: 4 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard       2.26 
## 2 mae     standard       1.62 
## 3 mape    standard      16.5  
## 4 rsq     standard       0.510

Ao calcularmos as métricas tanto no conjunto de treinamento quanto no de teste, podemos avaliar se o modelo está super-ajustando (overfitting):

reg_fit %>%
  augment(ring_train) %>%
  metrics(rings, .pred)
## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard       2.26 
## 2 rsq     standard       0.510
## 3 mae     standard       1.62
reg_fit %>%
  augment(ring_test) %>%
  metrics(rings, .pred)
## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard       2.24 
## 2 rsq     standard       0.505
## 3 mae     standard       1.63

7.2 Validação cruzada

Usaremos como exemplo 5 folds. Para fazer a reamostragem nos dados de treinamento, usaremos o comando vfold_cv:

set.seed(234)
ring_folds <- ring_train %>%
                vfold_cv(v = 5, strata = rings)
ring_folds
## #  5-fold cross-validation using stratification 
## # A tibble: 5 × 2
##   splits             id   
##   <list>             <chr>
## 1 <split [2670/670]> Fold1
## 2 <split [2672/668]> Fold2
## 3 <split [2672/668]> Fold3
## 4 <split [2673/667]> Fold4
## 5 <split [2673/667]> Fold5

Ajustando o modelo nas reamostras:

reg_cv <- reg_wf %>% fit_resamples(ring_folds)
reg_cv
## # Resampling results
## # 5-fold cross-validation using stratification 
## # A tibble: 5 × 4
##   splits             id    .metrics         .notes          
##   <list>             <chr> <list>           <list>          
## 1 <split [2670/670]> Fold1 <tibble [2 × 4]> <tibble [0 × 3]>
## 2 <split [2672/668]> Fold2 <tibble [2 × 4]> <tibble [0 × 3]>
## 3 <split [2672/668]> Fold3 <tibble [2 × 4]> <tibble [0 × 3]>
## 4 <split [2673/667]> Fold4 <tibble [2 × 4]> <tibble [0 × 3]>
## 5 <split [2673/667]> Fold5 <tibble [2 × 4]> <tibble [0 × 3]>

A função collect_metrics() extrai as métricas obtidas em cada reamostra e calcula a métrica da validação, geralmente através de uma média:

reg_cv %>%
  collect_metrics()
## # A tibble: 2 × 6
##   .metric .estimator  mean     n std_err .config             
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1 rmse    standard   2.33      5  0.126  Preprocessor1_Model1
## 2 rsq     standard   0.487     5  0.0323 Preprocessor1_Model1

Para calcular um conjunto escolhido de métricas, é preciso especificar o conjunto no argumento metrics dentro de fit_resamples:

reg_cv <- fit_resamples(reg_wf, 
                 ring_folds,
                 metrics = abalone_metrics)

reg_cv %>%
  collect_metrics()
## # A tibble: 4 × 6
##   .metric .estimator   mean     n std_err .config             
##   <chr>   <chr>       <dbl> <int>   <dbl> <chr>               
## 1 mae     standard    1.64      5  0.0377 Preprocessor1_Model1
## 2 mape    standard   16.7       5  0.295  Preprocessor1_Model1
## 3 rmse    standard    2.33      5  0.126  Preprocessor1_Model1
## 4 rsq     standard    0.487     5  0.0323 Preprocessor1_Model1

Através da validação cruzada, avaliamos o desempenho do modelo apenas com os dados de treinamento, sem usar os dados de teste.

A métrica obtida no conjunto de validação pode ser tomada como uma estimativa da métrica no conjunto de teste.

Caso seja necessário salvar as predições obtidas nas etapas de validação cruzada, para fazer um gráfico, por exemplo, usamos control_resamples:

ctrl_abalone <- control_resamples(save_pred = TRUE)
reg_cv <- fit_resamples(reg_wf, 
               ring_folds, 
               control = ctrl_abalone)
reg_preds <- collect_predictions(reg_cv)
reg_preds
## # A tibble: 3,340 × 5
##    .pred id     .row rings .config             
##    <dbl> <chr> <int> <dbl> <chr>               
##  1  8.03 Fold1     1     8 Preprocessor1_Model1
##  2  8.45 Fold1     3     8 Preprocessor1_Model1
##  3 10.9  Fold1    10     8 Preprocessor1_Model1
##  4  7.98 Fold1    23     7 Preprocessor1_Model1
##  5  7.26 Fold1    24     7 Preprocessor1_Model1
##  6  6.31 Fold1    25     7 Preprocessor1_Model1
##  7  5.86 Fold1    30     5 Preprocessor1_Model1
##  8  8.98 Fold1    34     8 Preprocessor1_Model1
##  9  9.51 Fold1    38     8 Preprocessor1_Model1
## 10  8.24 Fold1    39     7 Preprocessor1_Model1
## # ℹ 3,330 more rows

Podemos fazer um gráfico de preditos versus observados, separados por cada fold da validação cruzada:

reg_preds %>% 
  ggplot(aes(rings, .pred, color = id)) + 
  geom_abline(lty = 2, col = "gray", size = 1.5) +
  geom_point(alpha = 0.5) +
  coord_obs_pred()

7.3 Árvore de decisão

Um outro modelo a ser considerado é árvore de decisão. Vamos considerar o seguinte exemplo:

  tree_spec <- decision_tree(cost_complexity = 0.02) %>%
      set_mode("regression") %>%
      set_engine("rpart")

O modelo de árvore de regressão não requer pré-processamento, de forma que podemos usar um workflow com a fórmula e o modelo especificado, por exemplo:

tree_wf <- workflow(rings ~ ., tree_spec)
tree_wf
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: decision_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## rings ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Decision Tree Model Specification (regression)
## 
## Main Arguments:
##   cost_complexity = 0.02
## 
## Computational engine: rpart

E o ajuste é com os dados de treinamento:

tree_fit <- tree_wf %>% fit(data = ring_train)

Vamos então avaliar o desempenho da árvore de decisão usando validação cruzada:

tree_cv <- tree_wf %>% fit_resamples(ring_folds)
  
tree_cv %>% collect_metrics()
## # A tibble: 2 × 6
##   .metric .estimator  mean     n std_err .config             
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1 rmse    standard   2.50      5  0.0752 Preprocessor1_Model1
## 2 rsq     standard   0.402     5  0.0114 Preprocessor1_Model1

7.4 Conjunto de modelos

Quando queremos comparar vários modelos ao mesmo, fica muito trabalhoso fazer um por vez, como mostramos anteriormente.

Para isso, existe a função workflow_set() que gera um conjunto de workflows. Os argumentos desta função são:

  • preproc: formulas, receitas
  • models: modelos especificados usando parsnip
wf_set <- workflow_set(preproc = list(rec1 = rings ~ ., rec2 = ring_rec),
                       models = list(tree = tree_spec, reg = reg_spec),
                       cross = FALSE)

Agora podemos avaliar os modelos com as métricas desejadas usando a validação cruzada:

wf_set %>%
  workflow_map("fit_resamples", resamples = ring_folds) %>%
  rank_results()
## # A tibble: 4 × 9
##   wflow_id  .config         .metric  mean std_err     n preprocessor model  rank
##   <chr>     <chr>           <chr>   <dbl>   <dbl> <int> <chr>        <chr> <int>
## 1 rec2_reg  Preprocessor1_… rmse    2.33   0.126      5 recipe       line…     1
## 2 rec2_reg  Preprocessor1_… rsq     0.487  0.0323     5 recipe       line…     1
## 3 rec1_tree Preprocessor1_… rmse    2.50   0.0752     5 formula      deci…     2
## 4 rec1_tree Preprocessor1_… rsq     0.402  0.0114     5 formula      deci…     2

Se o argumento cross = TRUE o workflow_set faz um produto cruzado das receitas e modelos:

workflow_set(preproc = list(rec1 = rings ~ ., rec2 = ring_rec),
             models = list(tree = tree_spec, reg = reg_spec),
             cross = TRUE) %>%
  workflow_map("fit_resamples", resamples = ring_folds) %>%
  rank_results()
## # A tibble: 8 × 9
##   wflow_id  .config         .metric  mean std_err     n preprocessor model  rank
##   <chr>     <chr>           <chr>   <dbl>   <dbl> <int> <chr>        <chr> <int>
## 1 rec1_reg  Preprocessor1_… rmse    2.24   0.0921     5 formula      line…     1
## 2 rec1_reg  Preprocessor1_… rsq     0.521  0.0213     5 formula      line…     1
## 3 rec2_reg  Preprocessor1_… rmse    2.33   0.126      5 recipe       line…     2
## 4 rec2_reg  Preprocessor1_… rsq     0.487  0.0323     5 recipe       line…     2
## 5 rec1_tree Preprocessor1_… rmse    2.50   0.0752     5 formula      deci…     3
## 6 rec1_tree Preprocessor1_… rsq     0.402  0.0114     5 formula      deci…     3
## 7 rec2_tree Preprocessor1_… rmse    2.50   0.0752     5 recipe       deci…     4
## 8 rec2_tree Preprocessor1_… rsq     0.402  0.0114     5 recipe       deci…     4

Suponha que o modelo de regressão linear tenha sido o escolhido.

Vamos ajustar o modelo nos dados de treinamento e verificar o desempenho nos dados de teste.

Vimos os comandos fit() e predict()/augment(), mas podemos usar a função final_fit(), que combina esses passos.

final_fit <- last_fit(reg_wf, ring_split)

Lembre-se que o objeto ring_split tem informações sobre a separação dos dados originais em treino e teste.

final_fit
## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits             id               .metrics .notes   .predictions .workflow 
##   <list>             <chr>            <list>   <list>   <list>       <list>    
## 1 <split [3340/837]> train/test split <tibble> <tibble> <tibble>     <workflow>

Métricas calculadas para o conjunto de dados teste:

collect_metrics(final_fit) 
## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard       2.24  Preprocessor1_Model1
## 2 rsq     standard       0.505 Preprocessor1_Model1

Predições para o conjunto de dados teste:

collect_predictions(final_fit) %>%
  head()
## # A tibble: 6 × 5
##   .pred id                .row rings .config             
##   <dbl> <chr>            <int> <dbl> <chr>               
## 1  8.14 train/test split     2     7 Preprocessor1_Model1
## 2 10.6  train/test split     3     9 Preprocessor1_Model1
## 3  6.74 train/test split     5     7 Preprocessor1_Model1
## 4 11.2  train/test split    16    12 Preprocessor1_Model1
## 5  8.71 train/test split    19     7 Preprocessor1_Model1
## 6  8.69 train/test split    20     9 Preprocessor1_Model1
collect_predictions(final_fit) %>%
  ggplot(aes(rings, .pred)) + 
  geom_abline(lty = 2, col = "deeppink4", size = 1.5) +
  geom_point(alpha = 0.5) +
  coord_obs_pred()

Quais informações temos em final_fit?

extract_workflow(final_fit)
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 4 Recipe Steps
## 
## • step_normalize()
## • step_corr()
## • step_dummy()
## • step_poly()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:
## stats::lm(formula = ..y ~ ., data = data)
## 
## Coefficients:
##           (Intercept)                 height           shell_weight  
##              10.27917                0.59969                3.20871  
##            sex_infant               sex_male  shucked_weight_poly_1  
##              -1.02274               -0.02678             -127.72989  
## shucked_weight_poly_2  
##              -7.59663

Quais as estimativas dos parâmetros do modelo final?

final_fit %>% 
  extract_fit_parsnip() %>% 
  tidy()
## # A tibble: 7 × 5
##   term                   estimate std.error statistic   p.value
##   <chr>                     <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)             10.3       0.0746   138.    0        
## 2 height                   0.600     0.0704     8.52  2.35e- 17
## 3 shell_weight             3.21      0.0953    33.7   2.47e-214
## 4 sex_infant              -1.02      0.118     -8.65  7.83e- 18
## 5 sex_male                -0.0268    0.0964    -0.278 7.81e-  1
## 6 shucked_weight_poly_1 -128.        5.29     -24.2   6.75e-119
## 7 shucked_weight_poly_2   -7.60      2.50      -3.04  2.40e-  3

8 Hiperparâmetros

Algumas características dos modelos não podem ser estimadas diretamente dos dados.

Escolhemos um modelo de regressão linear, por exemplo, e usamos os dados de treinamento para obter os parâmetros do modelo.

No entanto, algumas escolhas são feitas antes do ajuste do modelo: usaremos alguma forma quadrática? interações? quais variáveis vamos considerar?

Algumas decisões devem ser feitas na etapa receita e outras devem ser feitas dentro do modelo.

Para ajuste fino, podemos testar workflows diferentes e avaliar o desempenho com validação cruzada.

8.1 Regressão linear com polinômio

Para um modelo de regressão linear em que uma das variáveis será considerada através de um polinômio, temos a seguinte receita:

ring_rec <-
  recipe(rings ~ ., data = ring_train) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_poly(shucked_weight, degree = tune())

Repare que, acima, não especificamos diretamente o grau do polinômio. Vamos escolher o melhor hiperparâmetro usando a função tune().

Com a receita definida, vamos agregar as informações:

regpol_wf <- workflow(ring_rec, linear_reg())
regpol_wf
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 2 Recipe Steps
## 
## • step_dummy()
## • step_poly()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear Regression Model Specification (regression)
## 
## Computational engine: lm

A função tune_grid() calcula um conjunto de métricas usando validação cruzada para avaliar o desempenho em um conjunto pré-determinado de hiperparâmetros de um modelo ou de uma receita:

set.seed(123)
regpol_res <- tune_grid(regpol_wf, ring_folds, grid = tibble(degree=1:6))
regpol_res
## # Tuning results
## # 5-fold cross-validation using stratification 
## # A tibble: 5 × 4
##   splits             id    .metrics          .notes          
##   <list>             <chr> <list>            <list>          
## 1 <split [2670/670]> Fold1 <tibble [12 × 5]> <tibble [0 × 3]>
## 2 <split [2672/668]> Fold2 <tibble [12 × 5]> <tibble [0 × 3]>
## 3 <split [2672/668]> Fold3 <tibble [12 × 5]> <tibble [0 × 3]>
## 4 <split [2673/667]> Fold4 <tibble [12 × 5]> <tibble [0 × 3]>
## 5 <split [2673/667]> Fold5 <tibble [12 × 5]> <tibble [0 × 3]>

Apresentando os resultados (média dos 5 folds) para cada grau de polinômio (hiperparâmetro) considerado:

collect_metrics(regpol_res)
## # A tibble: 12 × 7
##    degree .metric .estimator  mean     n std_err .config             
##     <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
##  1      1 rmse    standard   2.24      5  0.0921 Preprocessor1_Model1
##  2      1 rsq     standard   0.521     5  0.0213 Preprocessor1_Model1
##  3      2 rmse    standard   2.24      5  0.0952 Preprocessor2_Model1
##  4      2 rsq     standard   0.521     5  0.0222 Preprocessor2_Model1
##  5      3 rmse    standard   2.22      5  0.0866 Preprocessor3_Model1
##  6      3 rsq     standard   0.531     5  0.0191 Preprocessor3_Model1
##  7      4 rmse    standard   2.21      5  0.0845 Preprocessor4_Model1
##  8      4 rsq     standard   0.535     5  0.0187 Preprocessor4_Model1
##  9      5 rmse    standard   2.21      5  0.0813 Preprocessor5_Model1
## 10      5 rsq     standard   0.537     5  0.0176 Preprocessor5_Model1
## 11      6 rmse    standard   2.21      5  0.0811 Preprocessor6_Model1
## 12      6 rsq     standard   0.537     5  0.0176 Preprocessor6_Model1

Visualização gráfica usando autoplot():

autoplot(regpol_res, metric = "rmse")

Mostrando os melhores resultados:

show_best(regpol_res, metric = "rmse", n = 3)
## # A tibble: 3 × 7
##   degree .metric .estimator  mean     n std_err .config             
##    <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1      5 rmse    standard    2.21     5  0.0813 Preprocessor5_Model1
## 2      6 rmse    standard    2.21     5  0.0811 Preprocessor6_Model1
## 3      4 rmse    standard    2.21     5  0.0845 Preprocessor4_Model1

E se ajustarmos um outro modelo, também com algum outro hiperparâmetro? Como comparar?

8.2 Árvore de decisão

A seguir temos um workflow para árvore de decisão com ajuste de hiperparâmetro.

Primeiro, o modelo é especificado incluindo tune():

tree_spec <-
  decision_tree(
    cost_complexity = tune()
  ) %>%
  set_mode("regression") %>% 
  set_engine("rpart")

A receita é definida:

tree_rec <- 
  recipe(rings ~ ., data = ring_train) %>%
  step_dummy(all_nominal_predictors())

Essas duas informações são agregadaas com workflow:

tree_wf <- workflow(tree_rec, tree_spec)

Vamos usar tune_grid() para avaliar vários valores para o hiperparâmwtro:

set.seed(9)
tree_res <-
  tune_grid(tree_wf, 
            resamples = ring_folds, 
            grid = 15,
            metrics = abalone_metrics)

É possível fornecer um data.frame na opção grid, para ser mais específico.

Métricas obtidas através de validação cruzada considerando os valores de grid:

tree_res %>% collect_metrics()
## # A tibble: 60 × 7
##    cost_complexity .metric .estimator   mean     n std_err .config              
##              <dbl> <chr>   <chr>       <dbl> <int>   <dbl> <chr>                
##  1   0.00000000423 mae     standard    1.73      5  0.0492 Preprocessor1_Model01
##  2   0.00000000423 mape    standard   17.0       5  0.426  Preprocessor1_Model01
##  3   0.00000000423 rmse    standard    2.43      5  0.0858 Preprocessor1_Model01
##  4   0.00000000423 rsq     standard    0.462     5  0.0160 Preprocessor1_Model01
##  5   0.0000000452  mae     standard    1.73      5  0.0492 Preprocessor1_Model02
##  6   0.0000000452  mape    standard   17.0       5  0.426  Preprocessor1_Model02
##  7   0.0000000452  rmse    standard    2.43      5  0.0858 Preprocessor1_Model02
##  8   0.0000000452  rsq     standard    0.462     5  0.0160 Preprocessor1_Model02
##  9   0.000998      mae     standard    1.65      5  0.0502 Preprocessor1_Model03
## 10   0.000998      mape    standard   16.3       5  0.462  Preprocessor1_Model03
## # ℹ 50 more rows

Em resumo, o melhor resultado para regressão polinomial, segundo a validação cruzada:

regpol_res %>% 
  show_best(metric = "rmse", n = 1) %>% 
  select(.metric, .estimator, mean, n, std_err, .config)
## # A tibble: 1 × 6
##   .metric .estimator  mean     n std_err .config             
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1 rmse    standard    2.21     5  0.0813 Preprocessor5_Model1

E melhor resultado para regressão polinomial, segundo a validação cruzada:

tree_res %>% 
  show_best(metric = "rmse", n = 1) %>% 
  select(.metric, .estimator, mean, n, std_err, .config)
## # A tibble: 1 × 6
##   .metric .estimator  mean     n std_err .config              
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
## 1 rmse    standard    2.37     5  0.0941 Preprocessor1_Model15

Entre as duas opções (regressão polinomial e árvore), vamos ficar com a regressão polinomial. Vamos selecionar o modelo com hiperparâmtro de melhor desempenho:

best_rmse <- select_best(regpol_res, metric = "rmse")
best_rmse
## # A tibble: 1 × 2
##   degree .config             
##    <int> <chr>               
## 1      5 Preprocessor5_Model1

E ajustando o modelo final:

final_res <-
  regpol_wf %>% 
  finalize_workflow(best_rmse) %>%
  last_fit(ring_split)
final_res
## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits             id               .metrics .notes   .predictions .workflow 
##   <list>             <chr>            <list>   <list>   <list>       <list>    
## 1 <split [3340/837]> train/test split <tibble> <tibble> <tibble>     <workflow>

last_fit() ajusta o modelo final com os dados de treino e avalia o desempenho nos dados de teste.

Resultado no conjunto de teste:

final_res %>% collect_metrics()
## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard       2.18  Preprocessor1_Model1
## 2 rsq     standard       0.533 Preprocessor1_Model1

Workflow final

Guardar todos os passos do ajuste final (obtidos usando o conjunto de treinamento):

fitted_wf <- extract_workflow(final_res)
fitted_wf
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 2 Recipe Steps
## 
## • step_dummy()
## • step_poly()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:
## stats::lm(formula = ..y ~ ., data = data)
## 
## Coefficients:
##           (Intercept)                 length               diameter  
##              -0.05048               -6.66421                9.96106  
##                height           whole_weight         viscera_weight  
##               7.33184                9.41184               -8.93447  
##          shell_weight             sex_infant               sex_male  
##              10.52625               -0.96181               -0.01411  
## shucked_weight_poly_1  shucked_weight_poly_2  shucked_weight_poly_3  
##            -241.15746               -8.36663               24.57516  
## shucked_weight_poly_4  shucked_weight_poly_5  
##             -11.06719                8.32636

Obter valores preditos:

predict(fitted_wf, ring_test[1:3,])
## # A tibble: 3 × 1
##   .pred
##   <dbl>
## 1  8.06
## 2 11.4 
## 3  6.85

Workflow final:

final_res %>% 
  extract_fit_parsnip() %>%
  tidy()
## # A tibble: 14 × 5
##    term                   estimate std.error statistic  p.value
##    <chr>                     <dbl>     <dbl>     <dbl>    <dbl>
##  1 (Intercept)             -0.0505    0.740    -0.0682 9.46e- 1
##  2 length                  -6.66      2.23     -2.99   2.79e- 3
##  3 diameter                 9.96      2.59      3.84   1.24e- 4
##  4 height                   7.33      1.64      4.47   8.14e- 6
##  5 whole_weight             9.41      0.790    11.9    4.84e-32
##  6 viscera_weight          -8.93      1.42     -6.29   3.57e-10
##  7 shell_weight            10.5       1.25      8.44   4.57e-17
##  8 sex_infant              -0.962     0.115    -8.38   7.49e-17
##  9 sex_male                -0.0141    0.0926   -0.152  8.79e- 1
## 10 shucked_weight_poly_1 -241.       12.1     -19.9    3.33e-83
## 11 shucked_weight_poly_2   -8.37      3.79     -2.21   2.73e- 2
## 12 shucked_weight_poly_3   24.6       2.71      9.07   1.96e-19
## 13 shucked_weight_poly_4  -11.1       2.34     -4.73   2.34e- 6
## 14 shucked_weight_poly_5    8.33      2.25      3.70   2.23e- 4

Referências

Hastie, Trevor, Robert Tibshirani, and Jerome Friedman. 2001. The Elements of Statistical Learning. Springer Series in Statistics. New York, NY, USA: Springer New York Inc.
James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2013. An Introduction to Statistical Learning: With Applications in r. Springer. https://www.statlearning.com.
Kuhn, M., and J. Silge. 2022. Tidy Modeling with r. O’Reilly Media. https://www.tmwr.org.
Nash, Sellers, Warwick, and Wes Ford. 1995. “Abalone.” UCI Machine Learning Repository.