Uma Jornada pelos TidyModels em R - Classificação

Samara Kiihl1 and Tatiana Benaglia1

1Departamento de Estatística - IMECC - UNICAMP

← 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

Os dados estão disponíveis no Kaggle: https://www.kaggle.com/fmendes/diabetes-from-dat263x-lab01. Faça o download do arquivo diabetes.csv e salve na mesma pasta do arquivo .Rmd.

Vamos carregar o pacote tidyverse e ler o conjunto de dados diabetes:

library(tidyverse)
diabetes <- read_csv("diabetes.csv")

Algumas informações rápidas dos dados:

diabetes %>% glimpse()

## Rows: 15,000
## Columns: 10
## $ PatientID              <dbl> 1354778, 1147438, 1640031, 1883350, 1424119, 16…
## $ Pregnancies            <dbl> 0, 8, 7, 9, 1, 0, 0, 0, 8, 1, 1, 3, 5, 7, 0, 3,…
## $ PlasmaGlucose          <dbl> 171, 92, 115, 103, 85, 82, 133, 67, 80, 72, 88,…
## $ DiastolicBloodPressure <dbl> 80, 93, 47, 78, 59, 92, 47, 87, 95, 31, 86, 96,…
## $ TricepsThickness       <dbl> 34, 47, 52, 25, 27, 9, 19, 43, 33, 40, 11, 31, …
## $ SerumInsulin           <dbl> 23, 36, 35, 304, 35, 253, 227, 36, 24, 42, 58, …
## $ BMI                    <dbl> 43.50973, 21.24058, 41.51152, 29.58219, 42.6045…
## $ DiabetesPedigree       <dbl> 1.21319135, 0.15836498, 0.07901857, 1.28286985,…
## $ Age                    <dbl> 21, 23, 23, 43, 22, 26, 21, 26, 53, 26, 22, 23,…
## $ Diabetic               <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1,…

Temos \(n=15000\) observações.

Temos as seguintes variáveis:

• PatientID - identificador de cada paciente

• Pregnancies- número de gestações

• PlasmaGlucose - concentração de glicose no plasma sanguíneo medida duas horas após a ingestão de glicose

• DiastolicBloodPressure - pressão diastólica

• TricepsThickness - espessura da prega cutânea do tríceps

• SerumInsulin - nível de insulina circulante no soro medida duas horas após a ingestão de glicose.

• BMI - IMC

• DiabetesPedigree - representa uma pontuação de risco genético/familiar para diabetes mellitus (DM), baseada no histórico familiar do paciente.

• Age - idade em anos

• Diabetic - indicador de diabetes (1 se sim e 0 se não)

As 10 primeiras observações:

library(kableExtra)
diabetes %>% head(10) %>% kable(booktabs = TRUE)

PatientID	Pregnancies	PlasmaGlucose	DiastolicBloodPressure	TricepsThickness	SerumInsulin	BMI	DiabetesPedigree	Age	Diabetic
1354778	0	171	80	34	23	43.50973	1.2131914	21	0
1147438	8	92	93	47	36	21.24058	0.1583650	23	0
1640031	7	115	47	52	35	41.51152	0.0790186	23	0
1883350	9	103	78	25	304	29.58219	1.2828698	43	1
1424119	1	85	59	27	35	42.60454	0.5495419	22	0
1619297	0	82	92	9	253	19.72416	0.1034245	26	0
1660149	0	133	47	19	227	21.94136	0.1741598	21	0
1458769	0	67	87	43	36	18.27772	0.2361649	26	0
1201647	8	80	95	33	24	26.62493	0.4439474	53	1
1403912	1	72	31	40	42	36.88958	0.1039436	26	0

PatientID não será utilizada no modelo, portanto, vamos remover esta variável do conjunto de dados:

diabetes <- diabetes %>% select(-PatientID)

A variável resposta é Diabetic. Para problemas em que o objetivo é classificação, tidymodels espera que a variável resposta seja factor e considera o primeiro nível como o evento de interesse.

diabetes <- diabetes %>% 
  mutate(Diabetic = factor(Diabetic))

levels(diabetes$Diabetic)

## [1] "0" "1"

Aqui, queremos que o evento de interesse seja 1, que indica que o paciente é diabético, desta forma, usamos:

diabetes <- diabetes %>% 
  mutate(Diabetic = relevel(Diabetic, ref = "1"))

levels(diabetes$Diabetic)

## [1] "1" "0"

Algumas estatísticas sumárias:

summary(diabetes)

##   Pregnancies     PlasmaGlucose   DiastolicBloodPressure TricepsThickness
##  Min.   : 0.000   Min.   : 44.0   Min.   : 24.00         Min.   : 7.00   
##  1st Qu.: 0.000   1st Qu.: 84.0   1st Qu.: 58.00         1st Qu.:15.00   
##  Median : 2.000   Median :104.0   Median : 72.00         Median :31.00   
##  Mean   : 3.225   Mean   :107.9   Mean   : 71.22         Mean   :28.81   
##  3rd Qu.: 6.000   3rd Qu.:129.0   3rd Qu.: 85.00         3rd Qu.:41.00   
##  Max.   :14.000   Max.   :192.0   Max.   :117.00         Max.   :93.00   
##   SerumInsulin        BMI        DiabetesPedigree       Age        Diabetic 
##  Min.   : 14.0   Min.   :18.20   Min.   :0.07804   Min.   :21.00   1: 5000  
##  1st Qu.: 39.0   1st Qu.:21.26   1st Qu.:0.13774   1st Qu.:22.00   0:10000  
##  Median : 83.0   Median :31.77   Median :0.20030   Median :24.00            
##  Mean   :137.9   Mean   :31.51   Mean   :0.39897   Mean   :30.14            
##  3rd Qu.:195.0   3rd Qu.:39.26   3rd Qu.:0.61629   3rd Qu.:35.00            
##  Max.   :799.0   Max.   :56.03   Max.   :2.30159   Max.   :77.00

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)
tidymodels_prefer()

set.seed(1234)
diabetes_split <- initial_split(diabetes, prop = 0.8, strata = Diabetic)

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

diabetes_train <- training(diabetes_split)
diabetes_test <- testing(diabetes_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).

Podemos realizar análises exploratórias como de costume no conjunto de dados de treinamento. Por exemplo, qual é a distribuição da variável resposta para o conjunto de treinamento?

diabetes_train %>% 
  ggplot(aes(x = Diabetic)) +
  geom_bar()

A distribuição etária é semelhante entre as diabéticos e não diabéticos?

diabetes_train %>% 
  ggplot(aes(x = Age, fill = Diabetic, group = Diabetic)) +
  geom_density(position = "identity", alpha = .6)

A distribuição de PlasmaGlucose é semelhante entre as diabéticos e não diabéticos?

diabetes_train %>% 
  ggplot(aes(x = PlasmaGlucose, fill = Diabetic, group = Diabetic)) +
  geom_density(position = "identity", alpha = .6)

Podemos usar também boxplot para avaliar a distribuição de BMI entre diabéticos e não diabéticos:

diabetes_train %>%
  ggplot(aes(x = Diabetic, y = BMI, fill = Diabetic)) +
  geom_boxplot(alpha = 0.6) +
  labs(x = "Diabetes", y = "Índice de Massa Corporal (BMI)") +
  theme_minimal()

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, para especificar um modelo de regressão logística, usamos:

logistic_reg()

## Logistic Regression Model Specification (classification)
## 
## Computational engine: glm

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(logistic_reg)

## function (mode = "classification", engine = "glm", penalty = NULL, 
##     mixture = NULL) 
## NULL

Aqui vemos que o método default para logistic_reg() é glm (modelos lineares generalizados), mas outros modelos também estão disponíveis.

Se quisermos ajustar um modelo de regressão logística via modelos lineares generalizados regularizados por lasso e elastic-net, podemos especificar através da engine:

logistic_reg() %>% 
   set_engine("glmnet") 

## Logistic Regression Model Specification (classification)
## 
## Computational engine: glmnet

Todos os modelos disponíveis estã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:

rec_simple <- recipe(Diabetic ~ ., data = diabetes_train)

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

rec_simple %>% summary() %>% kable(booktabs = TRUE)

variable	type	role	source
Pregnancies	double , numeric	predictor	original
PlasmaGlucose	double , numeric	predictor	original
DiastolicBloodPressure	double , numeric	predictor	original
TricepsThickness	double , numeric	predictor	original
SerumInsulin	double , numeric	predictor	original
BMI	double , numeric	predictor	original
DiabetesPedigree	double , numeric	predictor	original
Age	double , numeric	predictor	original
Diabetic	factor , unordered, nominal	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:

diabetes_rec <- recipe(Diabetic ~ ., data = diabetes_train) %>%
  step_poly(Age, degree = 2)

A tabela a seguir apresenta os dados pré-processamos segundo a receita acima:

Diabetic	Pregnancies	PlasmaGlucose	DiastolicBloodPressure	TricepsThickness	SerumInsulin	BMI	DiabetesPedigree	Age_poly_1	Age_poly_2
0	0	171	80	34	23	43.510	1.213	-0.007	0.006
0	7	115	47	52	35	41.512	0.079	-0.005	0.002
0	0	82	92	9	253	19.724	0.103	-0.003	-0.003
0	0	133	47	19	227	21.941	0.174	-0.007	0.006
0	1	88	86	11	58	43.225	0.230	-0.006	0.004
0	3	94	96	31	36	21.294	0.259	-0.005	0.002
0	7	110	82	16	44	36.089	0.281	-0.004	-0.002
0	0	148	58	11	179	39.192	0.161	0.011	-0.012
0	0	92	84	8	324	21.866	0.258	0.002	-0.012
0	0	118	95	7	276	42.501	0.084	-0.005	0.000

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 logística:

logreg_spec <- logistic_reg() %>%
  set_engine("glm", family = "binomial") %>%
  set_mode("classification")

A receita a ser aplicada está no objeto diabetes_rec.

Vamos agregar essas duas informações no workflow:

logreg_wf <- workflow(diabetes_rec, logreg_spec)

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

logreg_fit <- logreg_wf %>% fit(data = diabetes_train)
logreg_fit

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_poly()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:  stats::glm(formula = ..y ~ ., family = ~"binomial", data = data)
## 
## Coefficients:
##            (Intercept)             Pregnancies           PlasmaGlucose  
##               6.928384               -0.278826               -0.008881  
## DiastolicBloodPressure        TricepsThickness            SerumInsulin  
##              -0.011326               -0.023160               -0.003905  
##                    BMI        DiabetesPedigree              Age_poly_1  
##              -0.052296               -0.987093              -79.422268  
##             Age_poly_2  
##              32.092675  
## 
## Degrees of Freedom: 11999 Total (i.e. Null);  11990 Residual
## Null Deviance:       15280 
## Residual Deviance: 10280     AIC: 10300

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

logreg_fit %>% tidy(conf.int = TRUE) %>% 
  kable(booktabs = TRUE)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	6.9283841	0.1918184	36.119500	0	6.5551673	7.3071513
Pregnancies	-0.2788264	0.0075937	-36.718150	0	-0.2937908	-0.2640212
PlasmaGlucose	-0.0088807	0.0007768	-11.432341	0	-0.0104063	-0.0073610
DiastolicBloodPressure	-0.0113256	0.0014974	-7.563761	0	-0.0142655	-0.0083954
TricepsThickness	-0.0231602	0.0017118	-13.529774	0	-0.0265251	-0.0198142
SerumInsulin	-0.0039049	0.0001836	-21.271690	0	-0.0042664	-0.0035467
BMI	-0.0522958	0.0026186	-19.971214	0	-0.0574471	-0.0471814
DiabetesPedigree	-0.9870933	0.0632462	-15.607145	0	-1.1113740	-0.8634180
Age_poly_1	-79.4222681	2.6000689	-30.546217	0	-84.5358127	-74.3425722
Age_poly_2	32.0926749	2.4744368	12.969689	0	27.2599637	36.9621785

logreg_fit %>% tidy(conf.int = TRUE, exponentiate = TRUE) %>% 
  kable(booktabs = TRUE)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	1.020843e+03	0.1918184	36.119500	0	7.028667e+02	1.490924e+03
Pregnancies	7.566713e-01	0.0075937	-36.718150	0	7.454324e-01	7.679573e-01
PlasmaGlucose	9.911586e-01	0.0007768	-11.432341	0	9.896476e-01	9.926660e-01
DiastolicBloodPressure	9.887383e-01	0.0014974	-7.563761	0	9.858358e-01	9.916398e-01
TricepsThickness	9.771059e-01	0.0017118	-13.529774	0	9.738236e-01	9.803808e-01
SerumInsulin	9.961027e-01	0.0001836	-21.271690	0	9.957427e-01	9.964596e-01
BMI	9.490481e-01	0.0026186	-19.971214	0	9.441719e-01	9.539143e-01
DiabetesPedigree	3.726583e-01	0.0632462	-15.607145	0	3.291065e-01	4.217182e-01
Age_poly_1	0.000000e+00	2.6000689	-30.546217	0	0.000000e+00	0.000000e+00
Age_poly_2	8.663066e+13	2.4744368	12.969689	0	6.900043e+11	1.128418e+16

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

logreg_fit %>% predict(diabetes_train)

## # A tibble: 12,000 × 1
##    .pred_class
##    <fct>      
##  1 0          
##  2 0          
##  3 0          
##  4 0          
##  5 0          
##  6 0          
##  7 0          
##  8 0          
##  9 0          
## 10 0          
## # ℹ 11,990 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, adicionando-os em uma coluna no conjunto de dados especificado:

logreg_fit %>% augment(diabetes_train)

## # A tibble: 12,000 × 12
##    .pred_class .pred_1 .pred_0 Pregnancies PlasmaGlucose DiastolicBloodPressure
##    <fct>         <dbl>   <dbl>       <dbl>         <dbl>                  <dbl>
##  1 0            0.289    0.711           0           171                     80
##  2 0            0.420    0.580           7           115                     47
##  3 0            0.0490   0.951           0            82                     92
##  4 0            0.0350   0.965           0           133                     47
##  5 0            0.0730   0.927           1            88                     86
##  6 0            0.0806   0.919           3            94                     96
##  7 0            0.351    0.649           7           110                     82
##  8 0            0.372    0.628           0           148                     58
##  9 0            0.146    0.854           0            92                     84
## 10 0            0.164    0.836           0           118                     95
## # ℹ 11,990 more rows
## # ℹ 6 more variables: TricepsThickness <dbl>, SerumInsulin <dbl>, BMI <dbl>,
## #   DiabetesPedigree <dbl>, Age <dbl>, Diabetic <fct>

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:

logreg_fit %>%
  augment(new_data = diabetes_train) %>% 
  head()

## # A tibble: 6 × 12
##   .pred_class .pred_1 .pred_0 Pregnancies PlasmaGlucose DiastolicBloodPressure
##   <fct>         <dbl>   <dbl>       <dbl>         <dbl>                  <dbl>
## 1 0            0.289    0.711           0           171                     80
## 2 0            0.420    0.580           7           115                     47
## 3 0            0.0490   0.951           0            82                     92
## 4 0            0.0350   0.965           0           133                     47
## 5 0            0.0730   0.927           1            88                     86
## 6 0            0.0806   0.919           3            94                     96
## # ℹ 6 more variables: TricepsThickness <dbl>, SerumInsulin <dbl>, BMI <dbl>,
## #   DiabetesPedigree <dbl>, Age <dbl>, Diabetic <fct>

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

• Matriz de Confusão
• Acurácia: \(\frac{TP+TN}{TP+FP+FN+TN}\)
• Sensibilidade: \(\frac{TP}{TP+FN}\)
• Especificidade: \(\frac{TN}{FP+TN}\)
• Brier score: métrica para modelos de classificação, análoga ao erro quadrático médio em modelos de regressão ${i=1}^n{k=1}^C(y_{ik}-_{ik})^2
• Kappa: métrica semelhante à acurácia, mas ajustada para o acaso. Ela compara a acurácia observada com a acurácia esperada caso as previsões fossem feitas aleatoriamente, sendo especialmente útil quando há desequilíbrio nas classes (ex.: uma classe muito mais frequente que as outras).

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

logreg_fit %>%
  augment(new_data = diabetes_train) %>%
  metrics(truth = Diabetic, estimate = .pred_class) %>% 
  kable(booktabs = TRUE)

.metric	.estimator	.estimate
accuracy	binary	0.7974167
kap	binary	0.5262132

Podemos especificar apenas uma métrica:

logreg_fit %>%
  augment(new_data = diabetes_train) %>%
  accuracy(truth = Diabetic, estimate = .pred_class) %>% 
  kable(booktabs = TRUE)

.metric	.estimator	.estimate
accuracy	binary	0.7974167

Matriz de confusão de treinamento:

logreg_fit %>%
  augment(new_data = diabetes_train) %>%
  conf_mat(truth = Diabetic, estimate = .pred_class)

##           Truth
## Prediction    1    0
##          1 2481  912
##          0 1519 7088

Um gráfico ilustrando os resultados:

logreg_fit %>%
  augment(new_data = diabetes_train) %>%
  conf_mat(truth = Diabetic, estimate = .pred_class) %>%
  autoplot(type = "heatmap")

Também é possível especificar um conjunto de métricas. No exemplo a seguir, diabetes_metrics é definido como um conjunto de três métricas de interesse:

diabetes_metrics <- metric_set(accuracy, sensitivity, specificity)

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

augment(logreg_fit, new_data = diabetes_train) %>%
  diabetes_metrics(truth = Diabetic, estimate = .pred_class) %>% 
  kable(booktabs = TRUE)

.metric	.estimator	.estimate
accuracy	binary	0.7974167
sensitivity	binary	0.6202500
specificity	binary	0.8860000

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

logreg_fit %>%
  augment(diabetes_train) %>%
  metrics(truth = Diabetic, estimate = .pred_class)

## # A tibble: 2 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.797
## 2 kap      binary         0.526

logreg_fit %>%
  augment(diabetes_test) %>%
  metrics(truth = Diabetic, estimate = .pred_class) 

## # A tibble: 2 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.784
## 2 kap      binary         0.492

Essa métrica usa como corte a probabilidade 0.5 para declarar um evento (se a probabilidade estimada for acima de 0.5, declaramos Diabético = 1 e 0 caso contrário).

E se decidirmos definir o corte como 0.8? A variação desse limite afeta as métricas de sensibilidade e especificidade.

Podemos usar a curva ROC (receiver operator characteristic):

• x-axis: taxa de falso positivo (1 - especificidade)
• y-axis: taxa de verdadeiro positivo (sensibilidade)

portanto, temos um cenário com sensibilidade e especificidade calculadas em todos os limites possíveis.

augment(logreg_fit, new_data = diabetes_train) %>% 
  roc_curve(truth = Diabetic, .pred_1) %>% 
    autoplot()

A area sob a curva ROC pode ser usada como métrica:

augment(logreg_fit, new_data = diabetes_train) %>% 
  roc_auc(truth = Diabetic, .pred_1) %>% 
  kable(booktabs = TRUE)

.metric	.estimator	.estimate
roc_auc	binary	0.8625254

e podemos fixar algumas linhas para visualizar os valores de corte (threshold) e seus respectivos TPR e FPR:

augment(logreg_fit, new_data = diabetes_train) %>% 
  roc_curve(truth = Diabetic, .pred_1) %>%
  slice(1, 500, 1000, 2000, 4000, 6000, 8000, 10000) %>% 
  kable(booktabs = TRUE)

.threshold	specificity	sensitivity
-Inf	0.000000	1.00000
0.0242410	0.062000	0.99950
0.0348151	0.124000	0.99850
0.0561580	0.244750	0.99000
0.1150272	0.474625	0.94975
0.2254710	0.678875	0.85825
0.4256154	0.846750	0.69400
0.6936874	0.959375	0.41925

7.2 Validação Cruzada

Usaremos 5 folds como exemplo. Para reamostrar os dados de treino, usaremos o comando vfold_cv:

set.seed(234)
diabetes_folds <- diabetes_train %>%
                vfold_cv(v = 5, strata = Diabetic)
diabetes_folds

## #  5-fold cross-validation using stratification 
## # A tibble: 5 × 2
##   splits              id   
##   <list>              <chr>
## 1 <split [9600/2400]> Fold1
## 2 <split [9600/2400]> Fold2
## 3 <split [9600/2400]> Fold3
## 4 <split [9600/2400]> Fold4
## 5 <split [9600/2400]> Fold5

Ajustando o modelo nas reamostragens:

logreg_cv <- logreg_wf %>% fit_resamples(diabetes_folds)
logreg_cv

## # Resampling results
## # 5-fold cross-validation using stratification 
## # A tibble: 5 × 4
##   splits              id    .metrics         .notes          
##   <list>              <chr> <list>           <list>          
## 1 <split [9600/2400]> Fold1 <tibble [3 × 4]> <tibble [0 × 3]>
## 2 <split [9600/2400]> Fold2 <tibble [3 × 4]> <tibble [0 × 3]>
## 3 <split [9600/2400]> Fold3 <tibble [3 × 4]> <tibble [0 × 3]>
## 4 <split [9600/2400]> Fold4 <tibble [3 × 4]> <tibble [0 × 3]>
## 5 <split [9600/2400]> Fold5 <tibble [3 × 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:

logreg_cv %>%
  collect_metrics() %>% 
  kable(booktabs = TRUE)

.metric	.estimator	mean	n	std_err	.config
accuracy	binary	0.7965833	5	0.0021707	Preprocessor1_Model1
brier_class	binary	0.1391407	5	0.0009101	Preprocessor1_Model1
roc_auc	binary	0.8618658	5	0.0015880	Preprocessor1_Model1

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

logreg_cv <- fit_resamples(logreg_wf, 
                 diabetes_folds,
                 metrics = diabetes_metrics)

logreg_cv %>%
  collect_metrics() %>% 
  kable(booktabs = TRUE)

.metric	.estimator	mean	n	std_err	.config
accuracy	binary	0.7965833	5	0.0021707	Preprocessor1_Model1
sensitivity	binary	0.6180000	5	0.0054715	Preprocessor1_Model1
specificity	binary	0.8858750	5	0.0025890	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_malaria <- control_resamples(save_pred = TRUE)

logreg_cv <- fit_resamples(logreg_wf, 
               diabetes_folds, 
               control = ctrl_malaria)

logreg_preds <- collect_predictions(logreg_cv)
logreg_preds

## # A tibble: 12,000 × 7
##    .pred_class .pred_1 .pred_0 id     .row Diabetic .config             
##    <fct>         <dbl>   <dbl> <chr> <int> <fct>    <chr>               
##  1 0            0.0735   0.927 Fold1     5 0        Preprocessor1_Model1
##  2 0            0.305    0.695 Fold1    20 0        Preprocessor1_Model1
##  3 0            0.109    0.891 Fold1    24 0        Preprocessor1_Model1
##  4 0            0.362    0.638 Fold1    37 0        Preprocessor1_Model1
##  5 0            0.0716   0.928 Fold1    41 0        Preprocessor1_Model1
##  6 0            0.148    0.852 Fold1    44 0        Preprocessor1_Model1
##  7 0            0.0791   0.921 Fold1    47 0        Preprocessor1_Model1
##  8 0            0.360    0.640 Fold1    49 0        Preprocessor1_Model1
##  9 0            0.0798   0.920 Fold1    50 0        Preprocessor1_Model1
## 10 0            0.280    0.720 Fold1    67 0        Preprocessor1_Model1
## # ℹ 11,990 more rows

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.005) %>%
    set_mode("classification") %>%
    set_engine("rpart", model = TRUE)
tree_spec

## Decision Tree Model Specification (classification)
## 
## Main Arguments:
##   cost_complexity = 0.005
## 
## Engine-Specific Arguments:
##   model = TRUE
## 
## Computational engine: rpart

O modelo de árvore de decisã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(Diabetic ~ ., tree_spec)
tree_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: decision_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Diabetic ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Decision Tree Model Specification (classification)
## 
## Main Arguments:
##   cost_complexity = 0.005
## 
## Engine-Specific Arguments:
##   model = TRUE
## 
## Computational engine: rpart

E ajustar com os dados de treinamento:

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

Visualização do modelo ajustado:

library(rpart.plot)

tree_fit %>% 
  extract_fit_engine() %>% 
  rpart.plot()

Vamos então avaliar o desempenho da árvore de decisão através da validação cruzada:

tree_cv <- tree_wf %>% fit_resamples(diabetes_folds)
  
tree_cv %>% collect_metrics() %>% 
  kable(booktabs = TRUE)

.metric	.estimator	mean	n	std_err	.config
accuracy	binary	0.9050833	5	0.0019913	Preprocessor1_Model1
brier_class	binary	0.0784368	5	0.0014673	Preprocessor1_Model1
roc_auc	binary	0.9335519	5	0.0022146	Preprocessor1_Model1

7.4 Floresta Aleatória

Podemos também considerar um modelo de floresta aleatória.

Aqui, usaremos set_engine e adicionaremos o argumento importance = "impurity", para que possamos ter pontuações de importância variáveis, para avaliarmos quais preditores são mais relevantes.

rf_spec <- rand_forest(trees = 1000) %>%
    set_mode("classification") %>% 
    set_engine("ranger", importance = "impurity")
rf_spec

## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   trees = 1000
## 
## Engine-Specific Arguments:
##   importance = impurity
## 
## Computational engine: ranger

Um workflow simples:

rf_wf <- workflow(Diabetic ~ ., rf_spec)
rf_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Diabetic ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   trees = 1000
## 
## Engine-Specific Arguments:
##   importance = impurity
## 
## Computational engine: ranger

Avaliação das métricas usando validação cruzada:

set.seed(2024) #RF uses random numbers, so we need to set the seed
rf_cv <- rf_wf %>% 
  fit_resamples(diabetes_folds)

rf_cv %>% collect_metrics() %>% 
  kable(booktabs = TRUE)

.metric	.estimator	mean	n	std_err	.config
accuracy	binary	0.9374167	5	0.0026333	Preprocessor1_Model1
brier_class	binary	0.0488073	5	0.0008318	Preprocessor1_Model1
roc_auc	binary	0.9824586	5	0.0007990	Preprocessor1_Model1

7.5 Conjunto de modelos

Quando queremos comparar vários modelos ao mesmo tempo, é muito trabalhoso fazer um de cada vez, como mostramos anteriormente.

tidymodels se mostra útil também nesta tarefa. Primeiro, podemos usar a função workflow_set() para gerar um conjunto de workflows que desejamos avaliar. Os argumentos desta função são:

• preproc: formulas, recipes
• models: models specified using parsnip

Aqui, definimos um conjunto de workflows de interesse:

wf_set <- workflow_set(preproc = list(rec1 = Diabetic ~ ., 
                                      rec2 = diabetes_rec, 
                                      rec1 = Diabetic ~ .),
                       models = list(tree = tree_spec, 
                                     logreg = logreg_spec, 
                                     rf = rf_spec),
                       cross = FALSE)

Para processar este conjunto de workflows, usamos a função workflow_map(). Podemos avaliar os modelos com as métricas desejadas usando a validação cruzada usando fit_resamples:

wf_set %>%
  workflow_map("fit_resamples", 
               resamples = diabetes_folds,
               metrics = diabetes_metrics,
               seed = 2024) %>%
  rank_results()

## # A tibble: 9 × 9
##   wflow_id    .config       .metric  mean std_err     n preprocessor model  rank
##   <chr>       <chr>         <chr>   <dbl>   <dbl> <int> <chr>        <chr> <int>
## 1 rec1_rf     Preprocessor… accura… 0.937 0.00263     5 formula      rand…     1
## 2 rec1_rf     Preprocessor… sensit… 0.892 0.00170     5 formula      rand…     1
## 3 rec1_rf     Preprocessor… specif… 0.960 0.00444     5 formula      rand…     1
## 4 rec1_tree   Preprocessor… accura… 0.905 0.00199     5 formula      deci…     2
## 5 rec1_tree   Preprocessor… sensit… 0.859 0.00805     5 formula      deci…     2
## 6 rec1_tree   Preprocessor… specif… 0.928 0.00447     5 formula      deci…     2
## 7 rec2_logreg Preprocessor… accura… 0.797 0.00217     5 recipe       logi…     3
## 8 rec2_logreg Preprocessor… sensit… 0.618 0.00547     5 recipe       logi…     3
## 9 rec2_logreg Preprocessor… specif… 0.886 0.00259     5 recipe       logi…     3

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

workflow_set(preproc = list(rec1 = Diabetic ~ ., 
                            rec2 = diabetes_rec),
             models = list(tree = tree_spec, 
                           logreg = logreg_spec, 
                           rf = rf_spec),
             cross = TRUE) %>%
  workflow_map("fit_resamples", 
               resamples = diabetes_folds,
               metrics = diabetes_metrics,
               seed = 2024) %>%
  rank_results()

## # A tibble: 18 × 9
##    wflow_id    .config      .metric  mean std_err     n preprocessor model  rank
##    <chr>       <chr>        <chr>   <dbl>   <dbl> <int> <chr>        <chr> <int>
##  1 rec2_rf     Preprocesso… accura… 0.938 0.00257     5 recipe       rand…     1
##  2 rec2_rf     Preprocesso… sensit… 0.897 0.00310     5 recipe       rand…     1
##  3 rec2_rf     Preprocesso… specif… 0.958 0.00444     5 recipe       rand…     1
##  4 rec1_rf     Preprocesso… accura… 0.937 0.00263     5 formula      rand…     2
##  5 rec1_rf     Preprocesso… sensit… 0.892 0.00170     5 formula      rand…     2
##  6 rec1_rf     Preprocesso… specif… 0.960 0.00444     5 formula      rand…     2
##  7 rec1_tree   Preprocesso… accura… 0.905 0.00199     5 formula      deci…     3
##  8 rec1_tree   Preprocesso… sensit… 0.859 0.00805     5 formula      deci…     3
##  9 rec1_tree   Preprocesso… specif… 0.928 0.00447     5 formula      deci…     3
## 10 rec2_tree   Preprocesso… accura… 0.905 0.00199     5 recipe       deci…     4
## 11 rec2_tree   Preprocesso… sensit… 0.859 0.00805     5 recipe       deci…     4
## 12 rec2_tree   Preprocesso… specif… 0.928 0.00447     5 recipe       deci…     4
## 13 rec2_logreg Preprocesso… accura… 0.797 0.00217     5 recipe       logi…     5
## 14 rec2_logreg Preprocesso… sensit… 0.618 0.00547     5 recipe       logi…     5
## 15 rec2_logreg Preprocesso… specif… 0.886 0.00259     5 recipe       logi…     5
## 16 rec1_logreg Preprocesso… accura… 0.789 0.00109     5 formula      logi…     6
## 17 rec1_logreg Preprocesso… sensit… 0.603 0.00617     5 formula      logi…     6
## 18 rec1_logreg Preprocesso… specif… 0.883 0.002       5 formula      logi…     6

Suponha que o modelo de regressão logística 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(logreg_wf, 
                      diabetes_split,
                      metrics = diabetes_metrics) 

Lembre-se que o objeto diabetes_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 [12000/3000]> train/test spl… <tibble> <tibble> <tibble>     <workflow>

Métricas calculadas para o conjunto de dados teste:

collect_metrics(final_fit) 

## # A tibble: 3 × 4
##   .metric     .estimator .estimate .config             
##   <chr>       <chr>          <dbl> <chr>               
## 1 accuracy    binary         0.784 Preprocessor1_Model1
## 2 sensitivity binary         0.59  Preprocessor1_Model1
## 3 specificity binary         0.880 Preprocessor1_Model1

Predições para o conjunto de dados teste:

collect_predictions(final_fit) %>%
  head()

## # A tibble: 6 × 5
##   .pred_class id                .row Diabetic .config             
##   <fct>       <chr>            <int> <fct>    <chr>               
## 1 0           train/test split     2 0        Preprocessor1_Model1
## 2 0           train/test split     5 0        Preprocessor1_Model1
## 3 0           train/test split     8 0        Preprocessor1_Model1
## 4 0           train/test split    10 0        Preprocessor1_Model1
## 5 1           train/test split    13 1        Preprocessor1_Model1
## 6 0           train/test split    19 0        Preprocessor1_Model1

Quais informações temos em final_fit?

extract_workflow(final_fit)

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_poly()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:  stats::glm(formula = ..y ~ ., family = ~"binomial", data = data)
## 
## Coefficients:
##            (Intercept)             Pregnancies           PlasmaGlucose  
##               6.928384               -0.278826               -0.008881  
## DiastolicBloodPressure        TricepsThickness            SerumInsulin  
##              -0.011326               -0.023160               -0.003905  
##                    BMI        DiabetesPedigree              Age_poly_1  
##              -0.052296               -0.987093              -79.422268  
##             Age_poly_2  
##              32.092675  
## 
## Degrees of Freedom: 11999 Total (i.e. Null);  11990 Residual
## Null Deviance:       15280 
## Residual Deviance: 10280     AIC: 10300

Quais as estimativas dos parâmetros do modelo final?

final_fit %>% 
  extract_fit_parsnip() %>% 
  tidy()

## # A tibble: 10 × 5
##    term                    estimate std.error statistic   p.value
##    <chr>                      <dbl>     <dbl>     <dbl>     <dbl>
##  1 (Intercept)              6.93     0.192        36.1  1.12e-285
##  2 Pregnancies             -0.279    0.00759     -36.7  3.75e-295
##  3 PlasmaGlucose           -0.00888  0.000777    -11.4  2.88e- 30
##  4 DiastolicBloodPressure  -0.0113   0.00150      -7.56 3.92e- 14
##  5 TricepsThickness        -0.0232   0.00171     -13.5  1.04e- 41
##  6 SerumInsulin            -0.00390  0.000184    -21.3  2.08e-100
##  7 BMI                     -0.0523   0.00262     -20.0  9.80e- 89
##  8 DiabetesPedigree        -0.987    0.0632      -15.6  6.51e- 55
##  9 Age_poly_1             -79.4      2.60        -30.5  6.35e-205
## 10 Age_poly_2              32.1      2.47         13.0  1.82e- 38

8 Hiperparâmetros (tuning)

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

Escolhemos um modelo de regressão logística, 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 logística com polinômio

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

diabetes_rec <-
  recipe(Diabetic ~ ., data = diabetes_train) %>%
  step_poly(Age, 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:

logregpol_wf <- workflow(diabetes_rec, logreg_spec)
logregpol_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_poly()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Logistic Regression Model Specification (classification)
## 
## Engine-Specific Arguments:
##   family = binomial
## 
## Computational engine: glm

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:

logregpol_res <- tune_grid(logregpol_wf, 
                           diabetes_folds, 
                           grid = tibble(degree=1:5))
logregpol_res

## # Tuning results
## # 5-fold cross-validation using stratification 
## # A tibble: 5 × 4
##   splits              id    .metrics          .notes          
##   <list>              <chr> <list>            <list>          
## 1 <split [9600/2400]> Fold1 <tibble [15 × 5]> <tibble [1 × 3]>
## 2 <split [9600/2400]> Fold2 <tibble [15 × 5]> <tibble [1 × 3]>
## 3 <split [9600/2400]> Fold3 <tibble [15 × 5]> <tibble [1 × 3]>
## 4 <split [9600/2400]> Fold4 <tibble [15 × 5]> <tibble [1 × 3]>
## 5 <split [9600/2400]> Fold5 <tibble [15 × 5]> <tibble [1 × 3]>
## 
## There were issues with some computations:
## 
##   - Warning(s) x5: glm.fit: fitted probabilities numerically 0 or 1 occurred
## 
## Run `show_notes(.Last.tune.result)` for more information.

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

collect_metrics(logregpol_res)

## # A tibble: 15 × 7
##    degree .metric     .estimator  mean     n  std_err .config             
##     <int> <chr>       <chr>      <dbl> <int>    <dbl> <chr>               
##  1      1 accuracy    binary     0.789     5 0.00109  Preprocessor1_Model1
##  2      1 brier_class binary     0.142     5 0.000820 Preprocessor1_Model1
##  3      1 roc_auc     binary     0.860     5 0.00144  Preprocessor1_Model1
##  4      2 accuracy    binary     0.797     5 0.00217  Preprocessor2_Model1
##  5      2 brier_class binary     0.139     5 0.000910 Preprocessor2_Model1
##  6      2 roc_auc     binary     0.862     5 0.00159  Preprocessor2_Model1
##  7      3 accuracy    binary     0.814     5 0.00254  Preprocessor3_Model1
##  8      3 brier_class binary     0.127     5 0.00168  Preprocessor3_Model1
##  9      3 roc_auc     binary     0.889     5 0.00244  Preprocessor3_Model1
## 10      4 accuracy    binary     0.818     5 0.00253  Preprocessor4_Model1
## 11      4 brier_class binary     0.124     5 0.00187  Preprocessor4_Model1
## 12      4 roc_auc     binary     0.894     5 0.00273  Preprocessor4_Model1
## 13      5 accuracy    binary     0.830     5 0.00285  Preprocessor5_Model1
## 14      5 brier_class binary     0.118     5 0.00138  Preprocessor5_Model1
## 15      5 roc_auc     binary     0.902     5 0.00208  Preprocessor5_Model1

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

autoplot(logregpol_res, metric = "accuracy")

Mostrando os melhores resultados:

show_best(logregpol_res, 
          metric = "accuracy", 
          n = 3)

## # A tibble: 3 × 7
##   degree .metric  .estimator  mean     n std_err .config             
##    <int> <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1      5 accuracy binary     0.830     5 0.00285 Preprocessor5_Model1
## 2      4 accuracy binary     0.818     5 0.00253 Preprocessor4_Model1
## 3      3 accuracy binary     0.814     5 0.00254 Preprocessor3_Model1

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

8.2 LASSO

Podemos ajustar uma regressão logística com regularização LASSO:

logreglasso_spec <- logistic_reg(penalty = tune(), mixture = 1) %>%
  set_engine("glmnet", family = "binomial") %>%
  set_mode("classification")

mixture = 1 especifica LASSO puro e o hiperparâmetro penalty será escolhido (ajuste fino) através da validação cruzada.

Não consideraremos o polinômio para Age neste workflow, então podemos configurá-lo diretamente com a fórmula e o modelo especificado:

logreglasso_wf <- workflow(Diabetic ~ ., logreglasso_spec)
logreglasso_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Diabetic ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = tune()
##   mixture = 1
## 
## Engine-Specific Arguments:
##   family = binomial
## 
## Computational engine: glmnet

Para avaliar o desempenho de múltiplos valores do hiperparâmetro usando validação cruzada:

set.seed(2024)
logreglasso_res <-
  tune_grid(logreglasso_wf, 
            resamples = diabetes_folds, 
            grid = 20,
            metrics = diabetes_metrics)

Métricas resultantes da validação cruzada considerando os valores do grid:

logreglasso_res %>% collect_metrics()

## # A tibble: 60 × 7
##     penalty .metric     .estimator  mean     n std_err .config              
##       <dbl> <chr>       <chr>      <dbl> <int>   <dbl> <chr>                
##  1 2.77e-10 accuracy    binary     0.789     5 0.00107 Preprocessor1_Model01
##  2 2.77e-10 sensitivity binary     0.6       5 0.00597 Preprocessor1_Model01
##  3 2.77e-10 specificity binary     0.883     5 0.00189 Preprocessor1_Model01
##  4 3.45e-10 accuracy    binary     0.789     5 0.00107 Preprocessor1_Model02
##  5 3.45e-10 sensitivity binary     0.6       5 0.00597 Preprocessor1_Model02
##  6 3.45e-10 specificity binary     0.883     5 0.00189 Preprocessor1_Model02
##  7 1.17e- 9 accuracy    binary     0.789     5 0.00107 Preprocessor1_Model03
##  8 1.17e- 9 sensitivity binary     0.6       5 0.00597 Preprocessor1_Model03
##  9 1.17e- 9 specificity binary     0.883     5 0.00189 Preprocessor1_Model03
## 10 5.42e- 9 accuracy    binary     0.789     5 0.00107 Preprocessor1_Model04
## # ℹ 50 more rows

autoplot(logreglasso_res)

O melhor resultado da validação cruzada, considerando a métrica de precisão:

show_best(logreglasso_res, metric = "accuracy", n = 1)

## # A tibble: 1 × 7
##    penalty .metric  .estimator  mean     n std_err .config              
##      <dbl> <chr>    <chr>      <dbl> <int>   <dbl> <chr>                
## 1 2.77e-10 accuracy binary     0.789     5 0.00107 Preprocessor1_Model01

8.3 Árvore de decisão

Abaixo está um workflow para uma árvore de decisão com ajuste de hiperparâmetro (cost_complexity). Primeiramente, o modelo é especificado incluindo tune():

tree_spec <-
  decision_tree(
    cost_complexity = tune()
  ) %>%
  set_mode("classification") %>% 
  set_engine("rpart", model = TRUE)

O workflow:

tree_wf <- workflow(Diabetic ~ ., tree_spec) 

Vamos usar tune_grid() para avaliar múltiplos valores para o hiperparâmetro:

tree_res <-
  tune_grid(tree_wf, 
            resamples = diabetes_folds, 
            grid = 30,
            metrics = diabetes_metrics)

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

Métricas obtidas por meio de validação cruzada considerando os valores da grade:

tree_res %>% collect_metrics()

## # A tibble: 90 × 7
##    cost_complexity .metric     .estimator  mean     n std_err .config           
##              <dbl> <chr>       <chr>      <dbl> <int>   <dbl> <chr>             
##  1        1.24e- 9 accuracy    binary     0.913     5 0.00159 Preprocessor1_Mod…
##  2        1.24e- 9 sensitivity binary     0.869     5 0.00393 Preprocessor1_Mod…
##  3        1.24e- 9 specificity binary     0.935     5 0.00306 Preprocessor1_Mod…
##  4        1.40e-10 accuracy    binary     0.913     5 0.00159 Preprocessor1_Mod…
##  5        1.40e-10 sensitivity binary     0.869     5 0.00393 Preprocessor1_Mod…
##  6        1.40e-10 specificity binary     0.935     5 0.00306 Preprocessor1_Mod…
##  7        2.22e- 4 accuracy    binary     0.914     5 0.00143 Preprocessor1_Mod…
##  8        2.22e- 4 sensitivity binary     0.874     5 0.00435 Preprocessor1_Mod…
##  9        2.22e- 4 specificity binary     0.935     5 0.00344 Preprocessor1_Mod…
## 10        5.45e- 9 accuracy    binary     0.913     5 0.00159 Preprocessor1_Mod…
## # ℹ 80 more rows

autoplot(tree_res)

O melhor resultado da validação cruzada, considerando a métrica de acurácia:

show_best(tree_res, metric = "accuracy", n = 1)

## # A tibble: 1 × 7
##   cost_complexity .metric  .estimator  mean     n std_err .config              
##             <dbl> <chr>    <chr>      <dbl> <int>   <dbl> <chr>                
## 1         0.00126 accuracy binary     0.920     5 0.00146 Preprocessor1_Model16

8.4 Floresta aleatória

Um hiperparâmetro comumente escolhido para modelos de floresta aleatória é mtry, que é o número de preditores que serão amostrados aleatoriamente em cada bifurcação ao criar os modelos de árvore.

Quando usamos o engine (algoritmo/método) ranger, o default é floor(sqrt(ncol(x))).

Vamos usar a validação cruzada para fazer o “ajuste fino” desse hiperparâmetro.

Começamos especificando o modelo e incluindo tune():

rf_spec <- rand_forest(trees = 1000,
                       mtry = tune()) %>%
    set_mode("classification") %>% 
    set_engine("ranger", importance = "impurity")
rf_spec

## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = 1000
## 
## Engine-Specific Arguments:
##   importance = impurity
## 
## Computational engine: ranger

O workflow:

rf_wf <- workflow(Diabetic ~ ., rf_spec)
rf_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Diabetic ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = 1000
## 
## Engine-Specific Arguments:
##   importance = impurity
## 
## Computational engine: ranger

Vamos usar tune_grid() para avaliar o desempenho de múltiplos valores do hiperparâmetro usando validação cruzada:

set.seed(2024)
rf_res <-
  tune_grid(rf_wf, 
            resamples = diabetes_folds, 
            grid = 15,
            metrics = diabetes_metrics)

Métricas resultantes da validação cruzada considerando os valores no grid:

rf_res %>% collect_metrics()

## # A tibble: 24 × 7
##     mtry .metric     .estimator  mean     n std_err .config             
##    <int> <chr>       <chr>      <dbl> <int>   <dbl> <chr>               
##  1     8 accuracy    binary     0.933     5 0.00207 Preprocessor1_Model1
##  2     8 sensitivity binary     0.894     5 0.00390 Preprocessor1_Model1
##  3     8 specificity binary     0.953     5 0.00432 Preprocessor1_Model1
##  4     3 accuracy    binary     0.938     5 0.00247 Preprocessor1_Model2
##  5     3 sensitivity binary     0.894     5 0.00258 Preprocessor1_Model2
##  6     3 specificity binary     0.960     5 0.00432 Preprocessor1_Model2
##  7     6 accuracy    binary     0.934     5 0.00262 Preprocessor1_Model3
##  8     6 sensitivity binary     0.892     5 0.00470 Preprocessor1_Model3
##  9     6 specificity binary     0.955     5 0.00437 Preprocessor1_Model3
## 10     7 accuracy    binary     0.934     5 0.00268 Preprocessor1_Model4
## # ℹ 14 more rows

autoplot(rf_res)

O melhor resultado da validação cruzada, considerando a acurácia como métrica:

show_best(rf_res, metric = "accuracy", n = 1)

## # A tibble: 1 × 7
##    mtry .metric  .estimator  mean     n std_err .config             
##   <int> <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1     3 accuracy binary     0.938     5 0.00247 Preprocessor1_Model2

8.5 Ajustando o modelo final (após tuning)

Vamos supor, por exemplo, que entre as quatro opções (regressão logística polinomial, LASSO, árvore de decisão e floresta aleatória), decidimos usar a regressão logística LASSO com melhor desempenho após a escolha do hiperparâmetro (tuning).

Precisamos então selecionar o modelo com o hiperparâmetro de melhor desempenho, usando a função select_best():

 best_acc <- select_best(logreglasso_res, metric = "accuracy")
 best_acc

## # A tibble: 1 × 2
##    penalty .config              
##      <dbl> <chr>                
## 1 2.77e-10 Preprocessor1_Model01

Para ajustar o modelo final, pegamos o workflow desejado (logregpol_wf neste exemplo) e usamos a função finalize_workflow() para especificar o hiperparâmetro com melhor desempenho. A função last_fit() ajusta este modelo final com os dados de treinamento e avalia o desempenho nos dados de teste.

final_lasso_fit <- logreglasso_wf %>% 
   finalize_workflow(best_acc) %>%
   last_fit(diabetes_split)
final_lasso_fit

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

Resultados no conjunto de teste:

final_lasso_fit %>% collect_metrics()

## # A tibble: 3 × 4
##   .metric     .estimator .estimate .config             
##   <chr>       <chr>          <dbl> <chr>               
## 1 accuracy    binary         0.776 Preprocessor1_Model1
## 2 roc_auc     binary         0.846 Preprocessor1_Model1
## 3 brier_class binary         0.150 Preprocessor1_Model1

Para o workflow final selecionado, podemos salvar todas as etapas de ajuste final (obtidas usando o conjunto de treinamento):

fitted_wf <- extract_workflow(final_lasso_fit)
fitted_wf

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Diabetic ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:  glmnet::glmnet(x = maybe_matrix(x), y = y, family = ~"binomial",      alpha = ~1) 
## 
##    Df  %Dev   Lambda
## 1   0  0.00 0.194600
## 2   1  2.25 0.177300
## 3   1  4.09 0.161500
## 4   2  6.46 0.147200
## 5   2  8.71 0.134100
## 6   2 10.59 0.122200
## 7   2 12.17 0.111300
## 8   3 13.70 0.101400
## 9   3 15.29 0.092430
## 10  4 16.91 0.084220
## 11  4 18.40 0.076730
## 12  4 19.68 0.069920
## 13  5 20.94 0.063710
## 14  5 22.09 0.058050
## 15  6 23.21 0.052890
## 16  6 24.22 0.048190
## 17  7 25.16 0.043910
## 18  7 26.02 0.040010
## 19  7 26.76 0.036460
## 20  7 27.39 0.033220
## 21  7 27.94 0.030270
## 22  8 28.43 0.027580
## 23  8 28.89 0.025130
## 24  8 29.29 0.022890
## 25  8 29.62 0.020860
## 26  8 29.91 0.019010
## 27  8 30.16 0.017320
## 28  8 30.38 0.015780
## 29  8 30.56 0.014380
## 30  8 30.71 0.013100
## 31  8 30.84 0.011940
## 32  8 30.95 0.010880
## 33  8 31.05 0.009911
## 34  8 31.13 0.009030
## 35  8 31.20 0.008228
## 36  8 31.25 0.007497
## 37  8 31.30 0.006831
## 38  8 31.34 0.006224
## 39  8 31.38 0.005671
## 40  8 31.40 0.005167
## 41  8 31.43 0.004708
## 42  8 31.45 0.004290
## 43  8 31.47 0.003909
## 44  8 31.48 0.003562
## 45  8 31.49 0.003245
## 46  8 31.50 0.002957
## 
## ...
## and 13 more lines.

Obtenha valores preditos para os dados de teste:

fitted_wf %>% augment(diabetes_test)

## # A tibble: 3,000 × 12
##    .pred_class .pred_1 .pred_0 Pregnancies PlasmaGlucose DiastolicBloodPressure
##    <fct>         <dbl>   <dbl>       <dbl>         <dbl>                  <dbl>
##  1 0            0.324    0.676           8            92                     93
##  2 0            0.107    0.893           1            85                     59
##  3 0            0.0389   0.961           0            67                     87
##  4 0            0.0624   0.938           1            72                     31
##  5 1            0.581    0.419           5           114                    101
##  6 0            0.288    0.712           8           117                     39
##  7 1            0.631    0.369           2            44                     81
##  8 0            0.136    0.864           0           119                     50
##  9 1            0.631    0.369           8           152                     83
## 10 0            0.170    0.830           2           118                     61
## # ℹ 2,990 more rows
## # ℹ 6 more variables: TricepsThickness <dbl>, SerumInsulin <dbl>, BMI <dbl>,
## #   DiabetesPedigree <dbl>, Age <dbl>, Diabetic <fct>

Resultados no conjunto de treinamento:

fitted_wf %>% 
  augment(diabetes_train) %>% 
  metrics(truth = Diabetic, estimate = .pred_class)

## # A tibble: 2 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.790
## 2 kap      binary         0.507

Para obter as estimativas do modelo:

final_lasso_fit %>%
  extract_fit_parsnip() %>%
  tidy() %>% 
  kable(booktabs = TRUE)

term	estimate	penalty
(Intercept)	8.4881010	0
Pregnancies	-0.2752126	0
PlasmaGlucose	-0.0084157	0
DiastolicBloodPressure	-0.0108160	0
TricepsThickness	-0.0224399	0
SerumInsulin	-0.0038390	0
BMI	-0.0514529	0
DiabetesPedigree	-0.9734048	0
Age	-0.0578455	0

8.6 Everything Everywhere All at Once

Podemos definir um conjunto de workflows e então avaliá-lo usando tune_grid.

wf_tune_set <- workflow_set(preproc = list(rec1 = Diabetic ~ ., 
                                           rec2 = diabetes_rec,
                                           rec1 = Diabetic ~ .,
                                           rec1 = Diabetic ~ .),
                            models = list(tree = tree_spec, 
                                          logreg = logreg_spec, 
                                          rf = rf_spec,
                                          lasso_log = logreglasso_spec),
                            cross = FALSE)

set.seed(2024)
wf_tune_res <-  wf_tune_set %>% 
  workflow_map(resamples = diabetes_folds, 
               grid = 40,
               metrics = diabetes_metrics,
               control = control_grid(save_pred = TRUE))

tune_results <- wf_tune_res %>% 
  collect_metrics()
tune_results

## # A tibble: 273 × 9
##    wflow_id  .config        preproc model .metric .estimator  mean     n std_err
##    <chr>     <chr>          <chr>   <chr> <chr>   <chr>      <dbl> <int>   <dbl>
##  1 rec1_tree Preprocessor1… formula deci… accura… binary     0.913     5 0.00148
##  2 rec1_tree Preprocessor1… formula deci… sensit… binary     0.869     5 0.00381
##  3 rec1_tree Preprocessor1… formula deci… specif… binary     0.935     5 0.00276
##  4 rec1_tree Preprocessor1… formula deci… accura… binary     0.913     5 0.00159
##  5 rec1_tree Preprocessor1… formula deci… sensit… binary     0.869     5 0.00393
##  6 rec1_tree Preprocessor1… formula deci… specif… binary     0.935     5 0.00306
##  7 rec1_tree Preprocessor1… formula deci… accura… binary     0.916     5 0.00202
##  8 rec1_tree Preprocessor1… formula deci… sensit… binary     0.857     5 0.00199
##  9 rec1_tree Preprocessor1… formula deci… specif… binary     0.945     5 0.00317
## 10 rec1_tree Preprocessor1… formula deci… accura… binary     0.887     5 0.00205
## # ℹ 263 more rows

Vamos selecionar o workflow com melhor desempenho, considerando a acurácia:

best_wf_id <- wf_tune_res %>% 
  rank_results() %>% 
  filter(.metric == "accuracy") %>%
  arrange(desc(mean)) %>%
  slice(1) %>%  # Select the best one
  pull(wflow_id)

best_wf_id

## [1] "rec1_rf"

Então, podemos extrair o melhor workflow da seguinte forma:

best_workflow <- wf_tune_set %>%
  extract_workflow(best_wf_id)

Para finalizar o workflow selecionado fixando o hiperparâmetro com melhor desempenho:

best_workflow_final <- best_workflow %>%
  finalize_workflow(
    wf_tune_res %>%
      extract_workflow_set_result(best_wf_id) %>%
      select_best(metric = "accuracy")
  )
best_workflow_final

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Diabetic ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = 2
##   trees = 1000
## 
## Engine-Specific Arguments:
##   importance = impurity
## 
## Computational engine: ranger

final_wf_fit <-
  best_workflow_final %>%
  last_fit(diabetes_split,
           metrics = diabetes_metrics)
final_wf_fit

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

Resultados no conjunto teste:

final_wf_fit %>% collect_metrics()

## # A tibble: 3 × 4
##   .metric     .estimator .estimate .config             
##   <chr>       <chr>          <dbl> <chr>               
## 1 accuracy    binary         0.937 Preprocessor1_Model1
## 2 sensitivity binary         0.894 Preprocessor1_Model1
## 3 specificity binary         0.958 Preprocessor1_Model1

Para o workflow final selecionado, salvamos todas as etapas de ajuste final (obtidas usando o conjunto de treinamento):

fitted_wf <- extract_workflow(final_wf_fit)
fitted_wf

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Diabetic ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Ranger result
## 
## Call:
##  ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~2L,      x), num.trees = ~1000, importance = ~"impurity", num.threads = 1,      verbose = FALSE, seed = sample.int(10^5, 1), probability = TRUE) 
## 
## Type:                             Probability estimation 
## Number of trees:                  1000 
## Sample size:                      12000 
## Number of independent variables:  8 
## Mtry:                             2 
## Target node size:                 10 
## Variable importance mode:         impurity 
## Splitrule:                        gini 
## OOB prediction error (Brier s.):  0.04772504

saveRDS(fitted_wf, "diabetes_rf_modelo.rds")

Valores preditos para o conjunto de teste:

final_wf_fit %>% augment()

## # A tibble: 3,000 × 10
##    .pred_class Pregnancies PlasmaGlucose DiastolicBloodPressure TricepsThickness
##    <fct>             <dbl>         <dbl>                  <dbl>            <dbl>
##  1 0                     8            92                     93               47
##  2 0                     1            85                     59               27
##  3 0                     0            67                     87               43
##  4 0                     1            72                     31               40
##  5 1                     5           114                    101               43
##  6 0                     8           117                     39               32
##  7 1                     2            44                     81               46
##  8 0                     0           119                     50               52
##  9 0                     8           152                     83               42
## 10 1                     2           118                     61               56
## # ℹ 2,990 more rows
## # ℹ 5 more variables: SerumInsulin <dbl>, BMI <dbl>, DiabetesPedigree <dbl>,
## #   Age <dbl>, Diabetic <fct>

Workflow final:

fitted_wf %>% 
  extract_fit_parsnip() 

## parsnip model object
## 
## Ranger result
## 
## Call:
##  ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~2L,      x), num.trees = ~1000, importance = ~"impurity", num.threads = 1,      verbose = FALSE, seed = sample.int(10^5, 1), probability = TRUE) 
## 
## Type:                             Probability estimation 
## Number of trees:                  1000 
## Sample size:                      12000 
## Number of independent variables:  8 
## Mtry:                             2 
## Target node size:                 10 
## Variable importance mode:         impurity 
## Splitrule:                        gini 
## OOB prediction error (Brier s.):  0.04772504

Visualização da importância dos preditores:

library(vip)
final_wf_fit %>% 
  extract_fit_parsnip() %>% 
  vip(num_features = 10)

9 Predição para novos dados

Com o modelo final escolhido e salvo, podemos utilizá-lo na predição de Diabetic quando tivermos acesso às variáveis preditoras de novos pacientes.

Para carregar o modelo salvo:

modelo_final <- readRDS("diabetes_rf_modelo.rds")

Preparando os dados de dois novos pacientes, por exemplo:

novos_pacientes <- tibble(Pregnancies = c(1, 2),
                        PlasmaGlucose = c(100, 180),
                        DiastolicBloodPressure = c(30, 80),
                        TricepsThickness = c( 15, 42),
                        SerumInsulin = c(39, 180),
                        BMI = c(25, 39),
                        DiabetesPedigree = c(0.15, 0.60),
                        Age = c(30, 50))

novos_pacientes

## # A tibble: 2 × 8
##   Pregnancies PlasmaGlucose DiastolicBloodPressure TricepsThickness SerumInsulin
##         <dbl>         <dbl>                  <dbl>            <dbl>        <dbl>
## 1           1           100                     30               15           39
## 2           2           180                     80               42          180
## # ℹ 3 more variables: BMI <dbl>, DiabetesPedigree <dbl>, Age <dbl>

e acrescentar as colunas com as predições (tanto de classe quanto de probabilidade para cada classe):

bind_cols(
  novos_pacientes,
  predict(modelo_final, new_data = novos_pacientes, type = "class"),
  predict(modelo_final, new_data = novos_pacientes, type = "prob")
)

## # A tibble: 2 × 11
##   Pregnancies PlasmaGlucose DiastolicBloodPressure TricepsThickness SerumInsulin
##         <dbl>         <dbl>                  <dbl>            <dbl>        <dbl>
## 1           1           100                     30               15           39
## 2           2           180                     80               42          180
## # ℹ 6 more variables: BMI <dbl>, DiabetesPedigree <dbl>, Age <dbl>,
## #   .pred_class <fct>, .pred_1 <dbl>, .pred_0 <dbl>

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.

Sharika, Tarin Sultana, Abdullah Al Farabe, Ghalib Ashraf, Nahian Raonak, and Amitabha Chakrabarty. 2021. “A Supervised Learning Approach by Machine Learning Algorithms to Predict Diabetes Mellitus (DM) Risk Score.” In Soft Computing: Theories and Applications, 289–300. Springer Singapore. https://doi.org/10.1007/978-981-16-1696-9_27.