Suport Vector Regression Machine - SVRM

Comparando com outros modelos

1 Carregando bibliotecas necessárias

library(tidymodels)

── Attaching packages ────────────────────────────────────── tidymodels 1.1.0 ──

✔ broom        1.0.5     ✔ recipes      1.0.6
✔ dials        1.2.0     ✔ rsample      1.1.1
✔ dplyr        1.1.2     ✔ tibble       3.2.1
✔ ggplot2      3.4.2     ✔ tidyr        1.3.0
✔ infer        1.0.4     ✔ tune         1.1.1
✔ modeldata    1.1.0     ✔ workflows    1.1.3
✔ parsnip      1.1.0     ✔ workflowsets 1.0.1
✔ purrr        1.0.1     ✔ yardstick    1.2.0

── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ purrr::discard() masks scales::discard()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ recipes::step()  masks stats::step()
• Use tidymodels_prefer() to resolve common conflicts.

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.4
✔ lubridate 1.9.2     ✔ stringr   1.5.0

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ readr::col_factor() masks scales::col_factor()
✖ purrr::discard()    masks scales::discard()
✖ dplyr::filter()     masks stats::filter()
✖ stringr::fixed()    masks recipes::fixed()
✖ dplyr::lag()        masks stats::lag()
✖ readr::spec()       masks yardstick::spec()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(GGally)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

library(skimr)

# Resolvendo possíveis conflitos entre o tidymodels e outras bibliotecas
tidymodels::tidymodels_prefer()

2 Importando a base de dados

Utilizando os dados de vinho vermelho🍷, disponíveis aqui, faça uma pequena análise exploratória dos dados. No link do Kaggle você consegue uma explicação sobre o que significa cada uma das variáveis.

Os dados, no meu caso, estão no diretório "../dados/winequality-red.csv". Você deverá alterar o path para o diretório encontra-se a base que deverá ser obtida no link acima.

dados <- readr::read_csv(file = "../dados/winequality-red.csv")

Rows: 1599 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (12): fixed acidity, volatile acidity, citric acid, residual sugar, chlo...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

3 Uma exploração rápida dos dados

É sempre importante olhar os dados antes de tentar modelar. Uma análise exploratória sempre será útil para identificarmos possíveis inconsistências.

visdat::vis_dat(dados)

O gráfico acima mostra que temos uma base de dados sem informações faltantes e todas as features presentes na base são numéricas. É uma situação confortável, haja vista que, aqui, não precisaremos nos preocupar com imputação de dados faltantes.

Um resumo dos dados poderá ser obtido utilizando a função glimpse do pacote dplyr que é carregado com a biblioteca tidyverse de R.

dados |> 
  dplyr::glimpse()

Rows: 1,599
Columns: 12
$ `fixed acidity`        <dbl> 7.4, 7.8, 7.8, 11.2, 7.4, 7.4, 7.9, 7.3, 7.8, 7…
$ `volatile acidity`     <dbl> 0.700, 0.880, 0.760, 0.280, 0.700, 0.660, 0.600…
$ `citric acid`          <dbl> 0.00, 0.00, 0.04, 0.56, 0.00, 0.00, 0.06, 0.00,…
$ `residual sugar`       <dbl> 1.9, 2.6, 2.3, 1.9, 1.9, 1.8, 1.6, 1.2, 2.0, 6.…
$ chlorides              <dbl> 0.076, 0.098, 0.092, 0.075, 0.076, 0.075, 0.069…
$ `free sulfur dioxide`  <dbl> 11, 25, 15, 17, 11, 13, 15, 15, 9, 17, 15, 17, …
$ `total sulfur dioxide` <dbl> 34, 67, 54, 60, 34, 40, 59, 21, 18, 102, 65, 10…
$ density                <dbl> 0.9978, 0.9968, 0.9970, 0.9980, 0.9978, 0.9978,…
$ pH                     <dbl> 3.51, 3.20, 3.26, 3.16, 3.51, 3.51, 3.30, 3.39,…
$ sulphates              <dbl> 0.56, 0.68, 0.65, 0.58, 0.56, 0.56, 0.46, 0.47,…
$ alcohol                <dbl> 9.4, 9.8, 9.8, 9.8, 9.4, 9.4, 9.4, 10.0, 9.5, 1…
$ quality                <dbl> 5, 5, 5, 6, 5, 5, 5, 7, 7, 5, 5, 5, 5, 5, 5, 5,…

É possível todas as correlações entre todas as variáveis da base, com a função data_vis_cor. Um gráfico útil com as correlações poderá ser obtido usando a função vis_cor, conforme abaixo:

visdat::data_vis_cor(dados)

# A tibble: 144 × 3
   row_1         row_2                  value
   <chr>         <chr>                  <dbl>
 1 fixed acidity fixed acidity         1     
 2 fixed acidity volatile acidity     -0.256 
 3 fixed acidity citric acid           0.672 
 4 fixed acidity residual sugar        0.115 
 5 fixed acidity chlorides             0.0937
 6 fixed acidity free sulfur dioxide  -0.154 
 7 fixed acidity total sulfur dioxide -0.113 
 8 fixed acidity density               0.668 
 9 fixed acidity pH                   -0.683 
10 fixed acidity sulphates             0.183 
# ℹ 134 more rows

visdat::vis_cor(dados)

Um gráfico de scatterplot para as variáveis numéricas poderá ser útil. Você poderá fazer isso, usando a função ggcatmat do pacote GGally.

dados |> 
  GGally::ggscatmat()

Warning: The dot-dot notation (`..scaled..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(scaled)` instead.
ℹ The deprecated feature was likely used in the GGally package.
  Please report the issue at <https://github.com/ggobi/ggally/issues>.

As bibliotecas GGally e skimr também possuem funções úteis que podem nos auxiliar no processo de exploração dos dados.

dados |> 
  GGally::ggpairs()

dados |> 
  skimr::skim()

Data summary

Name	dados
Number of rows	1599
Number of columns	12
_______________________
Column type frequency:
numeric	12
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
fixed acidity	0	1	8.32	1.74	4.60	7.10	7.90	9.20	15.90	▂▇▂▁▁
volatile acidity	0	1	0.53	0.18	0.12	0.39	0.52	0.64	1.58	▅▇▂▁▁
citric acid	0	1	0.27	0.19	0.00	0.09	0.26	0.42	1.00	▇▆▅▁▁
residual sugar	0	1	2.54	1.41	0.90	1.90	2.20	2.60	15.50	▇▁▁▁▁
chlorides	0	1	0.09	0.05	0.01	0.07	0.08	0.09	0.61	▇▁▁▁▁
free sulfur dioxide	0	1	15.87	10.46	1.00	7.00	14.00	21.00	72.00	▇▅▁▁▁
total sulfur dioxide	0	1	46.47	32.90	6.00	22.00	38.00	62.00	289.00	▇▂▁▁▁
density	0	1	1.00	0.00	0.99	1.00	1.00	1.00	1.00	▁▃▇▂▁
pH	0	1	3.31	0.15	2.74	3.21	3.31	3.40	4.01	▁▅▇▂▁
sulphates	0	1	0.66	0.17	0.33	0.55	0.62	0.73	2.00	▇▅▁▁▁
alcohol	0	1	10.42	1.07	8.40	9.50	10.20	11.10	14.90	▇▇▃▁▁
quality	0	1	5.64	0.81	3.00	5.00	6.00	6.00	8.00	▁▇▇▂▁

4 Construindo os workflows dos modelos

Iremos comparar os modelos de regressão linar utilizando elastic net, com o método $k$NN e suport vector regression machine - SVRM. Buscaremos pelo melhor modelo de cada uma das metodologias consideradas. Posteriormente iremos escolher o melhor modelo entre os melhores de cada uma das classes. A ideia é escolher o melhor modelo que consiga prever melhor a qualidade do vinho, i.e., prever a variável quality.

4.1 Divisão dos dados

Aqui usaremos as funções initial_split, training e testing para realizar o método hold-out (divisão inicial dos dados) em treino e teste. Vamos considerar $80\%$ para treino e $20\%$ para teste. A função initial_split irá realizar a divisão, porém, as funções training e testing são responsáveis para obtermos as tibbles da base de treino e teste, respectivamente.

Aqui, os dados está sendo estratificado pelo label, i.e., pela variável quality que desejamos prever:

set.seed(0) # Fixando uma semente
divisao_inicial <- rsample::initial_split(dados, prop = 0.8, strata = "quality")
treinamento <- rsample::training(divisao_inicial) # Conjunto de treinamento
teste <- rsample::testing(divisao_inicial) # Conjunto de teste

4.2 Tratamento dos dados (pré-processamento)

Apesar de não haver muito o que fazer nos dados que estamos utilizando nesse exemplo, em que nosso objetivo aqui é ter uma análise explicativa de como comparar modelos usando o tidymodels, iremos utilizar a biblioteca recipes. Os dados contém apenas variáveis numéricas com todas informações presentes, tornando o problema um pouco mais simples.

Na receita, iremos colocar as seguintes etapas:

1. Tomaremos o logarítmo de todas as variáveis peditoras (features);
2. Remover variáveis preditoras (features) que eventualmente estão altamente correlacionadas (usando a função step_corr);
3. Remover variáveis que possam ter variância próxima à zero, i.e., que sejam aproximadamente constantes (usando step_zv);
4. Normalizar os dados utilizando a função step_normalize.

Para que iremos remover variáveis altamente correlacionadas apenas nas variáveis preditoras, utilizamos a função all_predictors como argumentod a função step_corr. Já no passo de normalização dos dados, quando consideramos todas as variáveis numéricas, passamos para a função step_normalize a função all_numeric que especifica que deverá ser normalizado todas as variáveis numéricas. Na verdade, a normalização se dá em todas as variáveis numéricas, e portanto, esse argumento poderia ser omitido. Além disso, toda nossa base é formada por variáveis numéricas, o que torna redundante o uso, mas irei deixar explícito que todas as variáveis numéricas estão sendo normalizadas.

receita_1 <- 
  treinamento |> 
    recipe(formula = quality ~ .) |>
    step_YeoJohnson(all_predictors()) |>
    step_normalize(all_predictors()) |>
    step_zv(all_predictors()) |>
    step_corr(all_predictors())

receita_2 <- 
  treinamento |> 
    recipe(formula = quality ~ .) |>
    step_YeoJohnson(all_predictors()) |>
    step_normalize(all_predictors())

Como fazemos para observar se nosso pré-processamento funcionou?

Fácil, assim:

receita_1 |> 
  prep() |> 
  juice()

# A tibble: 1,278 × 12
   `fixed acidity` `volatile acidity` `citric acid` `residual sugar` chlorides
             <dbl>              <dbl>         <dbl>            <dbl>     <dbl>
 1         -0.446              1.00          -1.52            -0.602   -0.245 
 2         -0.162              1.77          -1.52             0.559    0.198 
 3         -0.162              1.28          -1.25             0.154    0.0774
 4         -0.446              1.00          -1.52            -0.602   -0.245 
 5         -0.446              0.813         -1.52            -0.845   -0.265 
 6         -0.0949             0.508         -1.12            -1.42    -0.386 
 7         -0.373             -0.0496         0.527            2.10    -0.346 
 8         -1.01               0.402         -0.994           -0.845    0.178 
 9         -0.373             -0.0496         0.527            2.10    -0.346 
10         -0.162              0.560          0.186           -1.42     0.521 
# ℹ 1,268 more rows
# ℹ 7 more variables: `free sulfur dioxide` <dbl>,
#   `total sulfur dioxide` <dbl>, density <dbl>, pH <dbl>, sulphates <dbl>,
#   alcohol <dbl>, quality <dbl>

A função prep estima uma receita de pré-processamento. Algumas funções step_* pode conter parâmetros que devem ser estimados. Inclusive, poderemos tunar esses parâmetros com a função tune do pacote tune. Por exemplo, se tivessemos interesse em interpolar uma feature usando a função step_ns, o argumento deg_free que refere-se ao grau do polinômio poderia ser “tunado”.

Todas as etapas de pré-processamento são estimadas em cima do conjunto de dados de treinamento. Por exemplo, na etapa em que realiza-se a normalização dos dados, a média a variância dos dados são estimadas uma única vez na base de dados completa, e sempre que essa refeita for aplicada a novos dados, será utilizado essa mesma média e variância, ou seja, não será recalculada com base no novo conjunto de dados.

Poderíamos utilizar a função bake ao invés da juice. A diferença de uma para a outra é que a função bake utiliza uma receita já estimada com prep e poderá ser aplicada à novos dados. Já a juice retorna a tibble com a receita preparada para o conjunto de dados de treinamento, ou ao conjunto de dados ao qual uma receita foi preparada com prep. Usando bake para o conjunto de dados de treinamento, poderíamos fazer:

receita_1 |> 
  prep() |> 
  bake(new_data = treinamento)

# A tibble: 1,278 × 12
   `fixed acidity` `volatile acidity` `citric acid` `residual sugar` chlorides
             <dbl>              <dbl>         <dbl>            <dbl>     <dbl>
 1         -0.446              1.00          -1.52            -0.602   -0.245 
 2         -0.162              1.77          -1.52             0.559    0.198 
 3         -0.162              1.28          -1.25             0.154    0.0774
 4         -0.446              1.00          -1.52            -0.602   -0.245 
 5         -0.446              0.813         -1.52            -0.845   -0.265 
 6         -0.0949             0.508         -1.12            -1.42    -0.386 
 7         -0.373             -0.0496         0.527            2.10    -0.346 
 8         -1.01               0.402         -0.994           -0.845    0.178 
 9         -0.373             -0.0496         0.527            2.10    -0.346 
10         -0.162              0.560          0.186           -1.42     0.521 
# ℹ 1,268 more rows
# ℹ 7 more variables: `free sulfur dioxide` <dbl>,
#   `total sulfur dioxide` <dbl>, density <dbl>, pH <dbl>, sulphates <dbl>,
#   alcohol <dbl>, quality <dbl>

4.3 Definindo os modelos

O código que segue faz a configuração realiza a configuração dos modelos que serão comparados. O código tune::tune() especifica que o respectivo parâmetro de sintonização será obtido no processo de validação cruzada, particularmente, um grid search.

modelo_elastic <- 
  parsnip::linear_reg(penalty = tune::tune(), mixture = tune::tune()) |> 
  parsnip::set_mode("regression") |> 
  parsnip::set_engine("glmnet")

modelo_knn <-
  parsnip::nearest_neighbor(
    neighbors = tune::tune(),
    dist_power = tune::tune(), 
    weight_func = "gaussian" 
  ) |> 
  parsnip::set_mode("regression") |> 
  parsnip::set_engine("kknn")

modelo_svm <- 
  parsnip::svm_rbf(
    cost = tune::tune(),
    rbf_sigma = tune::tune(),
    margin = tune::tune()
  ) |> 
  parsnip::set_mode("regression") |> 
  parsnip::set_engine("kernlab")

4.4 Criando o conjunto de validação

Uma vez que temos a divisão dos dados em conjunto de treinamento e conjunto de teste, precisamos construir o conjunto de validação, que necesse caso utilizaremos um $k$-fold cross-validation. Lembre-se que a validação cruzada ocorrerá sob a amostra de treinamento e necesso processo, a base de dados de treinamento será dividida em $k$ partes aproximadamente iguais em que trainamos o modelo em $k-1$ e validanos no fold restante, em que isso é feito $k$ vezes. Utilizando o conjunto $k-1$ em cada uma das divisões da validação cruzada em que diferentes combinações de hiperparâmetros são experimentadas e avaliada no conjunto de validação em cada divisão da validação cruzada.

O código que segue, em que utiliza a função vfold_cv da biblioteca rsample apenas cria a validação cruzada. Nesse caso, a validação será estratificada considerando a o label quality (qualidade do vinho), assim como foi feito na divisão inicial no hold-out (divisão inicial dos dados).

Utilizaremos $k=8$:

validacao_cruzada <- 
  treinamento |> 
  rsample::vfold_cv(v = 8L, strata = quality)

4.5 Criando um workflow completo

Aqui criaremos um workflow completo com todos modelos a serem comparados. Chamaremos ele de wf_todos:

wf_todos <-
  workflow_set(
    preproc = list(receita_1, receita_2),
    models = list(
      knn_fit = modelo_knn,
      elastic_fit = modelo_elastic,
      svm_fit = modelo_svm
    ),
    cross = TRUE
  )

Podemos manipular alguns argumentos que controlam aspectos da pesquisa em grade (grid search).. Fazemos isso com o uso da função control_grid da biblioteca parsnip. Por exemplo, podemos informar que desejamos paralelizar essa pesquisa, passando o argumento parallel_over = "resamples". Podemos também salvar as predições para cada um dos modelos especificando o argumento save_pred = TRUE. Se desejarmos anexar à saída o workflow, fazemos save_workflow = TRUE.

Assim, criaremos o objeto controle_grid, para que possamos passar a função workflow_map posteriormente. Tem-se:

controle_grid <- control_grid(
  save_pred = TRUE,
  save_workflow = TRUE,
  parallel_over = "resamples"
)

4.6 Trainando o modelo

Nessa etapa iremos treinar o modelo usando a função workflow_map do pacote workflow. Temos então:

treino <-
  wf_todos |> 
  workflow_map(
    resamples = validacao_cruzada,
    grid = 20L,
    control = controle_grid
  )

A função autoplot do pacote ggplot2 é útil para que possamos visualizar o desempenho de cada um dos modelos considerando a métrica do EQM. Isso é feito da seguinte forma:

autoplot(treino, metric = "rmse") + 
  labs(
    title = "Avaliação dos modelos de regressão",
    subtitle = "Utilizando a métrica do EQM"
  ) + 
  xlab("Rank dos Workflows") +
  ylab("Erro Quadrático Médio - EQM")

Perceba que no gráfico acima, temos os 8 modelos avaliados no $8$-folds cross-validation. Portanto, para cada um dos modelos comparados, temos 8 avaliações. Caso deseje avaliar as métricas do melhor cenário de cada um dos modelos, fazemos:

melhores <- 
  treino |> 
  rank_results(select_best = TRUE, rank_metric = "rmse")

autoplot(treino, select_best = TRUE)

Vamos agora selecionar o melhor modelo dentre os modelos comparados. Isso não quer dizer que o modelo seja bom para resolver o problema em questão. Para ver um rank e saber qual modelo e receita foram as melhores, fazemos:

treino |> 
  rank_results()

# A tibble: 240 × 9
   wflow_id         .config .metric  mean std_err     n preprocessor 
model  rank
   <chr>            <chr>   <chr>   <dbl>   
<dbl> <int> <chr>        <chr> <int>
 1 recipe_2_knn_fit Prepro… rmse    0.632  0.0168     8 recipe       
near…     1
 2 recipe_2_knn_fit Prepro… rsq     0.392  0.0182     8 recipe       
near…     1
 3 recipe_1_knn_fit Prepro… rmse    0.632  0.0168     8 recipe       
near…     2
 4 recipe_1_knn_fit Prepro… rsq     0.392  0.0182     8 recipe       
near…     2
 5 recipe_1_knn_fit Prepro… rmse    0.633  0.0165     8 recipe       
near…     3
 6 recipe_1_knn_fit Prepro… rsq     0.391  0.0174     8 recipe       
near…     3
 7 recipe_2_knn_fit Prepro… rmse    0.633  0.0165     8 recipe       
near…     4
 8 recipe_2_knn_fit Prepro… rsq     0.391  0.0174     8 recipe       
near…     4
 9 recipe_2_knn_fit Prepro… rmse    0.633  0.0169     8 recipe       
near…     5
10 recipe_2_knn_fit Prepro… rsq     0.390  0.0188     8 recipe       
near…     5
# ℹ 230 more rows

Assim, podemos perceber que a receita_1 combinada com o modelo $k$NN é o melhor escolha entre os modelos e receitas comparadas. Portanto:

melhor_modelo <- 
  treino |> 
  extract_workflow_set_result(id = "recipe_1_knn_fit") |> 
  select_best(metric = "rmse")

melhor_modelo

# A tibble: 1 × 3
  neighbors dist_power .config              
      <int>      <dbl> <chr>                
1        12       1.47 Preprocessor1_Model15

4.7 Avaliação final do melhor modelo

Após a escolha do melhor modelo e da estimação de seus hiperparâmetros, nesse caso, o modelo $k$NN com $k = 12$ e dist_power \approx 1.47, precisamos testar o desempenho do modelo final segundo a base de dados de teste. Para tanto, utilizamos a função last_fit do pacote tune. Temos que:

wf_final <- 
  treino |> 
  extract_workflow(id = "recipe_1_knn_fit") |> 
  finalize_workflow(melhor_modelo)

teste <- 
  wf_final |> 
  last_fit(split = divisao_inicial)

teste$.metrics

[[1]]
# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       0.629 Preprocessor1_Model1
2 rsq     standard       0.405 Preprocessor1_Model1

Agora que temos os hiperparâmetros estimados e temos uma boa estimativa do risco preditivo real do modelo final selecionado, poderemos preceder com um ajuste final, com toda a base de dados (treinamento + teste).

modelo_final <- 
  wf_final |> 
  fit(dados)

4.8 Salvando o modelo

Depois que temos o modelo finalizado, podemos ter alguns interesses de como utilizar o resultado, i.e., o modelo treinado. Entre alguns motivos, posso citar:

1. Salvar o modelo para uso no futuro, sem ter que retreinar;
2. Distribuir o modelo para que outras pessoas possam experimentar, sem terem que executar seu script R e retreinar o modelo;
3. Introduzir seu modelo treinado em uma API que irá consumir os resultados, i.e., consumir as previsões do modelo.

Nessas situações, é conveniente salvar o modelo em um arquivo serializado (R Data Serialization). Tais arquivos possuem a extensão .rds. Devemos fazer:

# Salvando o modelo em um arquivo serializado
saveRDS(modelo_final, file = "modelo_final.rds")

# Lendo o arquivo serializado com o modelo final
load(file = "modelo_final.rds")

# Fazendo novas previsões
predict(modelo_final, new_data = novos_dados)