Introdução Link para o cabeçalho

Neste tutorial iremos comparar o desempenho da regressão linear e das random forests na predição de preços de automóveis usados fabricados pela Volkswagen e disponíveis no Reino Unido. Para isso, utilizaremos o conjunto de dados 100,000 UK Used Car Data set. Todos os arquivos utilizados nesta análise podem ser baixados no github, inclusive os conjuntos de dados.

EDA Link para o cabeçalho

# pacotes necessarios

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::accumulate() masks foreach::accumulate()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ purrr::when()       masks foreach::when()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
theme_set(theme_bw())
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
# leitura dos dados

vw <- read_csv(file = "dados/vw.csv")
## Rows: 15157 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): model, transmission, fuelType
## dbl (6): year, price, mileage, tax, mpg, engineSize
## 
## ℹ 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.

Como em toda análise de dados, a primeira tarefa a ser realizada é a sua análise exploratória. Neste caso, quero verificar como está a relação entre as variáveis neste conjunto de dados, excetuando o modelo. Vou deixar esta variável de fora por enquanto, pois ela possui muitos níveis diferentes. Além disso, reordenei as colunas de modo que price fosse a primeira e, assim, a interpretação dos resultados ficasse mais simples. Note que, além disso, calculei a correlação de Spearman entre as variáveis, pois graficamente elas não apresentam linearidade ou normalidade.

vw %>%
  select(-model) %>%
  relocate(price) %>%
  relocate(where(is.character), .after = last_col()) %>%
  ggpairs(upper = list(continuous = wrap("cor", method = "spearman")))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Este gráfico nos dá uma ideia de como o preço está correlacionado com as outras variáveis preditoras. A maior correlação positiva do preço ocorre em relação ao ano de fabricação do veículo. Quanto mais recente o carro, maior o seu preço. O mesmo vale para o imposto (tax) e tamanho do motor (engineSize). Intuitivamente, nenhuma destas relações me parece estar longe do esperado.

O preço está correlacionado negativamente em relação à milhagem (mileage) e consumo de combustível (mpg, sigla para miles per galon). Novamente, de forma intuitiva, estes resultados prévios parecem estar condizentes com o esperado para dados assim.

A relação do preço com as variáveis qualitativas não está muito clara. A diferença de nível parece pequena e com muitos outliers. Entretanto, devido ao grande volume de dados (15157 observações), é possível que as diferenças de nível vistas nos boxplots sejam significativas.

Preparação dos Dados Link para o cabeçalho

Antes de partir para a modelagem em si, é necessário que os dados sejam preparados. Neste caso, a minha modelagem inicial vai contar com três etapas:

  1. Limpeza dos Dados
  2. Divisão dos Dados em Treinamento e Teste
  3. Pré-Processamento das Variáveis Preditoras

Limpeza dos Dados Link para o cabeçalho

Como é possível ver nos gráficos criados pela nossa análise exploratória, há uma quantidade maior de carros fabricados nos últimos 5 anos do que nos anos menos recentes. Abaixo veremos o quão grande é esta disparidade:

vw %>% 
  group_by(year) %>%
  count() %>%
  ungroup() %>%
  arrange(desc(year)) %>% 
  mutate(percentual = n/sum(n)*100, 
         acumulado = cumsum(percentual)) %>%
  print(n = Inf)
## # A tibble: 21 × 4
##     year     n percentual acumulado
##    <dbl> <int>      <dbl>     <dbl>
##  1  2020  1046    6.90         6.90
##  2  2019  4669   30.8         37.7 
##  3  2018  1509    9.96        47.7 
##  4  2017  2947   19.4         67.1 
##  5  2016  2647   17.5         84.6 
##  6  2015  1153    7.61        92.2 
##  7  2014   580    3.83        96.0 
##  8  2013   315    2.08        98.1 
##  9  2012    80    0.528       98.6 
## 10  2011    57    0.376       99.0 
## 11  2010    41    0.271       99.3 
## 12  2009    31    0.205       99.5 
## 13  2008    27    0.178       99.6 
## 14  2007    20    0.132       99.8 
## 15  2006    16    0.106       99.9 
## 16  2005     8    0.0528      99.9 
## 17  2004     3    0.0198      99.9 
## 18  2003     2    0.0132     100.  
## 19  2002     1    0.00660    100.  
## 20  2001     4    0.0264     100.  
## 21  2000     1    0.00660    100

Como é possível ver, 92,2% dos carros foram fabricados de 2015 em diante. Ou seja, manter carros mais antigos do que esta data no conjunto de dados vai causar um desbalanceamento muito grande, prejudicandos os resultados. Além disso, se criarmos uma tabela que conte o número de modelos de carro por ano, a situação fica bem complicada:

vw %>% 
  select(year, model) %>%
  group_by(year, model) %>%
  count() %>%
  ungroup() %>%
  complete(year, model, fill = list(n = 0)) %>%
  pivot_wider(names_from = "year", values_from = "n")
## # A tibble: 27 × 22
##    model   `2000` `2001` `2002` `2003` `2004` `2005` `2006` `2007` `2008` `2009`
##    <chr>    <int>  <int>  <int>  <int>  <int>  <int>  <int>  <int>  <int>  <int>
##  1 Amarok       0      0      0      0      0      0      0      0      0      0
##  2 Arteon       0      0      0      0      0      0      0      0      0      0
##  3 Beetle       0      3      0      1      0      1      3      3      0      1
##  4 Caddy        0      0      0      0      0      0      0      0      0      0
##  5 Caddy …      0      0      0      0      0      0      0      0      0      0
##  6 Caddy …      0      0      0      0      0      0      0      0      0      0
##  7 Caddy …      0      0      0      0      0      0      0      0      0      0
##  8 Califo…      0      0      0      0      0      0      0      0      0      0
##  9 Carave…      0      0      0      0      0      1      2      0      0      0
## 10 CC           0      0      0      0      0      0      0      0      0      0
## # ℹ 17 more rows
## # ℹ 11 more variables: `2010` <int>, `2011` <int>, `2012` <int>, `2013` <int>,
## #   `2014` <int>, `2015` <int>, `2016` <int>, `2017` <int>, `2018` <int>,
## #   `2019` <int>, `2020` <int>

Há uma quantidade de zeros muito grande na tabela construída. Graficamente, a situação é bastante clara:

vw %>% 
  select(year, model) %>%
  group_by(year, model) %>%
  count() %>%
  ungroup() %>%
  complete(year, model, fill = list(n = 0)) %>%
  mutate(faltante = ifelse(n == 0, "Sim", "Não")) %>%
  ggplot(., aes(x = year, y = model, fill = faltante)) +
  geom_tile(colour = "white") +
  labs(x = "Ano de Fabricação", 
       y = "Modelo do Carro", 
       fill = "Dado faltante?") +
  scale_fill_viridis_d()

A grande quantidade de retângulos amarelos no gráfico acima mostra como há muitos dados faltantes. Em particular, 59,79% das combinações entre modelo de carro e ano não possuem nenhuma ocorrência. Por isso, baseado no fato de que 92,2% dos carros foram fabricados de 2015 em diante, vou restringir arbitrariamente este conjunto de dados para apenas os anos mais recentes, a partir desta data. Com isso, temos o seguinte:

vw_bk <- vw

vw <- 
  vw_bk %>%
  filter(year >= 2015)

vw %>% 
  select(year, model) %>%
  group_by(year, model) %>%
  count() %>%
  ungroup() %>%
  complete(year, model, fill = list(n = 0)) %>%
  mutate(faltante = ifelse(n == 0, "Sim", "Não")) %>%
  ggplot(., aes(x = year, y = model, fill = faltante)) +
  geom_tile(colour = "white") +
  scale_x_continuous(breaks = seq(2015, 2020), labels = seq(2015, 2020)) +
  labs(x = "Ano de Fabricação", 
       y = "Modelo do Carro", 
       fill = "Dado faltante?") +
  scale_fill_viridis_d()

Com isso, apenas 25,64% das combinações estão faltando. Ainda assim, modelos como Eos, Caddy Maxi e Caddy Life parecem ter um comportamento não ideal, ainda com dados faltantes para muitos anos. Por isso, vou retirar todas as ocorrências deles do banco de dados.

vw <- 
  vw %>%
  filter(model != "Eos" & model != "Caddy Maxi" & model != "Caddy Life")

vw %>% 
  select(year, model) %>%
  group_by(year, model) %>%
  count() %>%
  ungroup() %>%
  complete(year, model, fill = list(n = 0)) %>%
  mutate(faltante = ifelse(n == 0, "Sim", "Não")) %>%
  ggplot(., aes(x = year, y = model, fill = faltante)) +
  geom_tile(colour = "white") +
  scale_x_continuous(breaks = seq(2015, 2020), labels = seq(2015, 2020)) +
  labs(x = "Ano de Fabricação", 
       y = "Modelo do Carro", 
       fill = "Dado faltante?") +
  scale_fill_viridis_d()

Agora sim, com apenas 20,29% de combinações entre anos e modelos de carros faltantes, podemos proceder com a análise. Além disso, continuamos com 13960 linhas no conjunto de dados a ser analisado, uma perda de apenas 7.9% em relação às 15157 linhas do conjunto de dados original.

Divisão dos Dados em Treinamento e Teste Link para o cabeçalho

Um post antigo deste blog, intitulado Tutorial: Como Fazer o Seu Primeiro Projeto de Data Science, discute a razão do porquê dividirmos os dados em treino e teste. Por este motivo, não irei justificar novamente a razão de termos que fazer algo nesta linha.

A partir de agora, utilizaremos o pacote tidymodels para tratar os dados antes das análises preditivas. Vamos dividir os dados em 75% para o treinamento e 25% para o teste.

# pacote necessario

library(tidymodels)
## ── Attaching packages ────────────────────────────────────── tidymodels 1.1.1 ──
## ✔ broom        1.0.5      ✔ rsample      1.2.0 
## ✔ dials        1.2.1      ✔ tune         1.1.2 
## ✔ infer        1.0.6      ✔ workflows    1.1.4 
## ✔ modeldata    1.3.0      ✔ workflowsets 1.0.1 
## ✔ parsnip      1.2.0      ✔ yardstick    1.3.0 
## ✔ recipes      1.0.10
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ purrr::accumulate() masks foreach::accumulate()
## ✖ scales::discard()   masks purrr::discard()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ recipes::fixed()    masks stringr::fixed()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ yardstick::spec()   masks readr::spec()
## ✖ recipes::step()     masks stats::step()
## ✖ purrr::when()       masks foreach::when()
## • Use suppressPackageStartupMessages() to eliminate package startup messages
# 75% dos dados como treino

set.seed(2410)

vw_split <- initial_split(vw, prop = .75)

# criar os conjuntos de dados de treino e teste

vw_treino <- training(vw_split)
vw_teste  <- testing(vw_split)

Pré-Processamento das Variáveis Preditoras Link para o cabeçalho

Com os dados limpos e separados em treinamento e teste, o próximo passo para a análise é pré-processar as variáveis. São dois os passos que serão executados a partir de agora. O primeiro é a criação de variáveis dummy. Estas variáveis permitem que dados qualitativos sejam analisados em um contexto de regressão. Basicamente transformamos dados categóricos em numéricos.

O segundo passo é centrar e escalar as variáveis preditoras. Isto é, deixá-las todas com média 0 e desvio padrão 1. Este passo ajuda na convergência dos algoritmos que utilizaremos para nossa modelagem.

Ao final, teremos um objeto chamado vw_rec, que lembrará quais são os passos que utilizamos para preparar os dados para análise, de modo que as mesmas transformações possam ser aplicadas tanto no conjunto de treino quanto no conjunto de teste.

vw_rec <- 
  # regressao que serah aplicada aos dados
  recipe(price ~ ., 
         data = vw_treino) %>%
  # criacao das variaveis dummy
  step_dummy(where(is.character)) %>%
  # center/scale
  step_center(where(is.numeric), -price) %>% 
  step_scale(where(is.numeric), -price) %>% 
  # funcao para aplicar a transformacao aos dados
  prep()

Em seguida, a receita vw_rec é aplicada aos conjuntos de treino e teste.

# aplicar a transformacao aos dados

vw_treino_t <- juice(vw_rec)

# preparar o conjunto de teste

vw_teste_t <- bake(vw_rec,
                   new_data = vw_teste)

Com os dados preparados, é hora de partirmos para a modelagem propriamente dita.

Modelagem Link para o cabeçalho

O objetivo a partir de agora, com os dados preparados e transformados, é prever o preço dos carros em função das outras variáveis presentes no conjunto de dados. Duas abordagens serão utilizadas e comparadas. A primeira delas é a regressão linear. Em seguida, vou utilizar random forests. Para que ambas sejam comparáveis, vou utilizar o pacote tidymodels, de modo que as métricas sejam as mesmas para ambas as abordagens.

Regressão Linear Link para o cabeçalho

Embora tenhamos visto que os dados não possuem relações lineares entre si, vou começar com a regressão linear para ter uma baseline com a qual trabalhar.

O pacote tidymodels exige que definamos uma série de objetos intermediários para que o modelo possa ser ajutado aos dados. Mais informações sobre o significado de cada um deles podem ser encontradas em meu curso de Introdução à Modelagem de Big Data ou no próprio site do pacote tidymodels. Como o bjetivo deste texto é ser apenas uma aplicação das técnicas deste pacote, deixarei estas explicações para lá.

Na ordem, iremos criar o modelo a ser ajustado, o workflow para o ajuste, a definição da validação cruzada e ajustaremos o modelo definido com estas características.

# modelo

vw_lm <- 
  linear_reg() %>% 
  set_engine("lm")

vw_lm
## Linear Regression Model Specification (regression)
## 
## Computational engine: lm
# criar workflow

vw_wflow <- 
  workflow() %>% 
  add_model(vw_lm) %>% 
  add_formula(price ~ .)

# definicao da validacao cruzada

set.seed(2410)

vw_treino_cv <- vfold_cv(vw_treino, v = 10)

vw_treino_cv
## #  10-fold cross-validation 
## # A tibble: 10 × 2
##    splits              id    
##    <list>              <chr> 
##  1 <split [9423/1047]> Fold01
##  2 <split [9423/1047]> Fold02
##  3 <split [9423/1047]> Fold03
##  4 <split [9423/1047]> Fold04
##  5 <split [9423/1047]> Fold05
##  6 <split [9423/1047]> Fold06
##  7 <split [9423/1047]> Fold07
##  8 <split [9423/1047]> Fold08
##  9 <split [9423/1047]> Fold09
## 10 <split [9423/1047]> Fold10
# modelo ajustado com validacao cruzada

vw_lm_fit_cv <- fit_resamples(vw_wflow, vw_treino_cv)

vw_lm_fit_cv
## # Resampling results
## # 10-fold cross-validation 
## # A tibble: 10 × 4
##    splits              id     .metrics         .notes          
##    <list>              <chr>  <list>           <list>          
##  1 <split [9423/1047]> Fold01 <tibble [2 × 4]> <tibble [0 × 3]>
##  2 <split [9423/1047]> Fold02 <tibble [2 × 4]> <tibble [0 × 3]>
##  3 <split [9423/1047]> Fold03 <tibble [2 × 4]> <tibble [0 × 3]>
##  4 <split [9423/1047]> Fold04 <tibble [2 × 4]> <tibble [0 × 3]>
##  5 <split [9423/1047]> Fold05 <tibble [2 × 4]> <tibble [0 × 3]>
##  6 <split [9423/1047]> Fold06 <tibble [2 × 4]> <tibble [0 × 3]>
##  7 <split [9423/1047]> Fold07 <tibble [2 × 4]> <tibble [0 × 3]>
##  8 <split [9423/1047]> Fold08 <tibble [2 × 4]> <tibble [0 × 3]>
##  9 <split [9423/1047]> Fold09 <tibble [2 × 4]> <tibble [0 × 3]>
## 10 <split [9423/1047]> Fold10 <tibble [2 × 4]> <tibble [0 × 3]>
# resultados

collect_metrics(vw_lm_fit_cv)
## # A tibble: 2 × 6
##   .metric .estimator     mean     n  std_err .config             
##   <chr>   <chr>         <dbl> <int>    <dbl> <chr>               
## 1 rmse    standard   2505.       10 30.7     Preprocessor1_Model1
## 2 rsq     standard      0.893    10  0.00289 Preprocessor1_Model1

Como as métricas da validação cruzada estão satisfatórias, iremos ajustar o modelo no conjunto de treino completo para, em seguida, verificar o resultado obtido no conjunto de teste.

# ajuste do modelo no conjunto de treino completo

vw_lm_fit_treino <- fit(vw_wflow, vw_treino)

# resultados no conjunto de teste

resultado <- 
  vw_teste %>%
  bind_cols(predict(vw_lm_fit_treino, vw_teste) %>%
              rename(predicao_lm = .pred))

# resultado final

metrics(resultado, 
        truth = price, 
        estimate = predicao_lm)
## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard    2413.   
## 2 rsq     standard       0.893
## 3 mae     standard    1785.

Comparando as métricas obtidas no conjunto de treino com o conjunto de teste, é possível dizer que não há sobreajuste, dado que os resultados estão muito próximos.

Ao comparar os valores preditos pelo modelo com os que estavam no conjunto teste, obtemos o gráfico a seguir.

# grafico final

ggplot(resultado, aes(x = price, y = predicao_lm)) +
  geom_point() +
  labs(x = "Valores Observados", y = "Valores Preditos") +
  geom_abline(intercept = 0, slope = 1) +
  coord_fixed()

Caso a nossa predição tivesse ficado perfeita, o gráfico acima deveria se comportar como uma reta. Como é possível perceber, o comportamento não está como o ideal. Por isso, tentaremos uma nova abordagem para os dados. Iremos utilizar um algoritmo muito famoso para isso, chamado random forests.

Random Forests Link para o cabeçalho

Uma grande vantagem das random forests sobre a regressão linear é a inexistência de uam exigência a respeito da linearidade dos dados. Com isso, torna-se algo muito mais versátil do que um método que só pode lidar com dados lineares.

Como estamos utilizando o pacote tidymodels, é possível reaproveitar grande parte do código utilizado para ajustar a regressão linear. As funções a serem utilizadas são rigorosamente as mesmas, com pequenas mudanças em seus argumentos.

Entretanto, será necessário adicionar algumas funções extras, a fim de realizar o tuning dos parâmetros para o modelo ajustado, a fim de minimizar o erro de predição do preço dos automóveis.

# definicao do tuning

vw_rf_tune <-
  rand_forest(
    mtry = tune(),
    trees = 1000,
    min_n = tune()
  ) %>%
  set_mode("regression") %>%
  set_engine("ranger", importance = "impurity")

# grid de procura

vw_rf_grid <- grid_regular(mtry(range(1, 8)),
                           min_n(range(10, 50)))

vw_rf_grid
## # A tibble: 9 × 2
##    mtry min_n
##   <int> <int>
## 1     1    10
## 2     4    10
## 3     8    10
## 4     1    30
## 5     4    30
## 6     8    30
## 7     1    50
## 8     4    50
## 9     8    50
# workflow

vw_rf_tune_wflow <- 
  workflow() %>%
  add_model(vw_rf_tune) %>%
  add_formula(price ~ .)

# definicao da validacao cruzada

set.seed(2410)

vw_treino_cv <- vfold_cv(vw_treino_t, v = 10)

# avaliacao do modelo

vw_rf_fit_tune <- 
  vw_rf_tune_wflow %>% 
  tune_grid(
    resamples = vw_treino_cv,
    grid = vw_rf_grid
  )

# resultados

collect_metrics(vw_rf_fit_tune)
## # A tibble: 18 × 8
##     mtry min_n .metric .estimator     mean     n  std_err .config             
##    <int> <int> <chr>   <chr>         <dbl> <int>    <dbl> <chr>               
##  1     1    10 rmse    standard   2408.       10 41.1     Preprocessor1_Model1
##  2     1    10 rsq     standard      0.921    10  0.00232 Preprocessor1_Model1
##  3     4    10 rmse    standard   1645.       10 34.3     Preprocessor1_Model2
##  4     4    10 rsq     standard      0.954    10  0.00183 Preprocessor1_Model2
##  5     8    10 rmse    standard   1703.       10 32.1     Preprocessor1_Model3
##  6     8    10 rsq     standard      0.950    10  0.00175 Preprocessor1_Model3
##  7     1    30 rmse    standard   2448.       10 39.9     Preprocessor1_Model4
##  8     1    30 rsq     standard      0.918    10  0.00239 Preprocessor1_Model4
##  9     4    30 rmse    standard   1692.       10 34.6     Preprocessor1_Model5
## 10     4    30 rsq     standard      0.951    10  0.00188 Preprocessor1_Model5
## 11     8    30 rmse    standard   1723.       10 32.1     Preprocessor1_Model6
## 12     8    30 rsq     standard      0.949    10  0.00174 Preprocessor1_Model6
## 13     1    50 rmse    standard   2503.       10 47.6     Preprocessor1_Model7
## 14     1    50 rsq     standard      0.914    10  0.00273 Preprocessor1_Model7
## 15     4    50 rmse    standard   1753.       10 36.0     Preprocessor1_Model8
## 16     4    50 rsq     standard      0.948    10  0.00195 Preprocessor1_Model8
## 17     8    50 rmse    standard   1781.       10 33.7     Preprocessor1_Model9
## 18     8    50 rsq     standard      0.946    10  0.00183 Preprocessor1_Model9

Novamente, as métricas da validação cruzada estão satisfatórias. Agora iremos verificar o desempenho dos hiperparâmetros do modelo, a fim de verificar se foram encontrados pontos de mínimo para o RMSE em cada um deles.

vw_rf_fit_tune %>%
  collect_metrics() %>%
  mutate(min_n = factor(min_n)) %>%
  ggplot(., aes(x = mtry, y = mean, colour = min_n, group = min_n)) +
  geom_line() +
  geom_point() +
  facet_grid(.metric ~ ., scales = "free") +
  scale_x_continuous(breaks = seq(1, 9, 2)) +
  scale_colour_viridis_d()

Tudo certo de acordo com o gráfico acima. Conseguimos obter o mínimo de RMSE (ou seja, o menor erro possível) e o máximo de \(R^2\) (nosso pseudo coeficiente de determinação) de maneira global. Os resultados apresentados no gráfico poderiam ser refinados se fosse criado um objeto vw_rf_tune com mais pontos para mtry, min_n ou até mesmo trees. Entretanto, isso faria com que mais modelos fossem gerados, o que elevaria o tempo de execução do ajuste dos modelos considerados.

Como foram gerados vários modelos diferentes a partir dos hiperparâmetros definidos, é preciso agora encontrar aqueles modelos que geraram os melhores resultados. Após, iremos ajustar no conjunto de treino completo aquele modelo cujos hiperparâmetros geraram o menor RMSE.

# melhores modelos

vw_rf_fit_tune %>%
  show_best("rmse")
## # A tibble: 5 × 8
##    mtry min_n .metric .estimator  mean     n std_err .config             
##   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1     4    10 rmse    standard   1645.    10    34.3 Preprocessor1_Model2
## 2     4    30 rmse    standard   1692.    10    34.6 Preprocessor1_Model5
## 3     8    10 rmse    standard   1703.    10    32.1 Preprocessor1_Model3
## 4     8    30 rmse    standard   1723.    10    32.1 Preprocessor1_Model6
## 5     4    50 rmse    standard   1753.    10    36.0 Preprocessor1_Model8
# melhor modelo

vw_rf_best <- 
  vw_rf_fit_tune %>%
  select_best("rmse")

vw_rf_final <-
  vw_rf_tune_wflow %>%
  finalize_workflow(vw_rf_best)

vw_rf_final <- fit(vw_rf_final, 
                   vw_treino_t)

vw_rf_final
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## price ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Ranger result
## 
## Call:
##  ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~4L,      x), num.trees = ~1000, min.node.size = min_rows(~10L, x),      importance = ~"impurity", num.threads = 1, verbose = FALSE,      seed = sample.int(10^5, 1)) 
## 
## Type:                             Regression 
## Number of trees:                  1000 
## Sample size:                      10470 
## Number of independent variables:  8 
## Mtry:                             4 
## Target node size:                 10 
## Variable importance mode:         impurity 
## Splitrule:                        variance 
## OOB prediction error (MSE):       2658260 
## R squared (OOB):                  0.9545377

Agora que o melhor modelo foi encontrado, podemos verificar os resultados obtidos por ele no conjunto de teste.

# resultados no conjunto de teste

resultado_rf <- 
  vw_teste_t %>%
  bind_cols(predict(vw_rf_final, vw_teste_t) %>%
              rename(predicao_rf = .pred))

metrics(resultado_rf, 
        truth = price, 
        estimate = predicao_rf)
## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard    1513.   
## 2 rsq     standard       0.958
## 3 mae     standard    1006.
metrics(resultado, 
        truth = price, 
        estimate = predicao_lm)
## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard    2413.   
## 2 rsq     standard       0.893
## 3 mae     standard    1785.

Como podemos ver, o RMSE no conjunto de teste para o random forest ficou em 1514 (1630 no conjunto de treino, indicando que não houve sobreajuste), bem abaixo dos 2413 obtidos pela regressão linear. O mesmo vale para o \(R^2\) do random forest, bem acima daquele obtido pela regressão linear.

Ao comparar graficamente os valores preditos pelo modelo com aqueles disponíveis no conjunto de teste, temos um resultado muito bom e, como era de se supor, bastante superior ao da regressão linear. A variância destes dados está muito menor, indicando que o ajuste foi muito bem feito.

# grafico final

ggplot(resultado_rf, aes(x = price, y = predicao_rf)) +
  geom_point() +
  labs(x = "Valores Observados", y = "Valores Preditos") +
  geom_abline(intercept = 0, slope = 1) +
  coord_fixed()

Por fim, podemos encontrar as variáveis mais importantes para a determinação do preço dos automóveis.

# importancia das variaveis

library(vip)
## 
## Attaching package: 'vip'
## The following object is masked from 'package:utils':
## 
##     vi
vw_rf_final %>% 
  extract_fit_parsnip() %>% 
  vip(scale = TRUE)

A variável que melhor determina o preço dos automóveis fabricados pela Volkswagen no mercado de usados no Reino Unido é o consumo, seguido de longe pelo tamanho do motor, ano de fabricação, modelo do carro e assim por diante. O tipo de câmbio é o que menos importa neste mercado.

Conclusões Link para o cabeçalho

Utilizar o tidymodels para ajustar modelos aos nossos dados é bastante simples. A grande vantagem é poder usar um conjunto padrão de funções que podem ser reaproveitadas pelos mais diferentes métodos de modelagem.

Fica a cargo do leitor ajustar modelos para os outros fabricantes de automóveis disponíveis neste conjunto de dados.

Todos os arquivos utilizados podem ser baixados no github.