13  Poszukiwanie optymalnego modelu

W przypadku projektów z nowymi zbiorami danych, które nie zostały jeszcze dobrze poznane, osoba zajmująca się danymi może być zmuszona do sprawdzenia wielu kombinacji modeli i różnych kombinacji metod przygotowania danych. Powszechne jest to, że nie posiadamy wiedzy a priori na temat tego, która metoda będzie działać najlepiej z nowym zestawem danych.

Dobrą strategią jest poświęcenie trochę uwagi na wypróbowanie różnych podejść do modelowania, określenie, co działa najlepiej, a następnie zainwestowanie dodatkowego czasu na dostosowanie / optymalizację małego zestawu modeli.

Aby zademonstrować, jak przesiewać zestaw wielu modeli, użyjemy jako przykładu danych mieszanki betonowej z książki Applied Predictive Modeling (Khun i Johnson 2013). W rozdziale 10 tej książki zademonstrowano modele do przewidywania wytrzymałości na ściskanie mieszanek betonowych z wykorzystaniem składników jako zmiennych niezależnych. Oceniono wiele różnych modeli z różnymi zestawami predyktorów i typami metod wstępnego przetwarzania.

13.1 Przygotowanie danych

Kod
library(tidymodels)
library(finetune)
tidymodels_prefer()
data(concrete, package = "modeldata")
glimpse(concrete)
Rows: 1,030
Columns: 9
$ cement               <dbl> 540.0, 540.0, 332.5, 332.5, 198.6, 266.0, 380.0, …
$ blast_furnace_slag   <dbl> 0.0, 0.0, 142.5, 142.5, 132.4, 114.0, 95.0, 95.0,…
$ fly_ash              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ water                <dbl> 162, 162, 228, 228, 192, 228, 228, 228, 228, 228,…
$ superplasticizer     <dbl> 2.5, 2.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,…
$ coarse_aggregate     <dbl> 1040.0, 1055.0, 932.0, 932.0, 978.4, 932.0, 932.0…
$ fine_aggregate       <dbl> 676.0, 676.0, 594.0, 594.0, 825.5, 670.0, 594.0, …
$ age                  <int> 28, 28, 270, 365, 360, 90, 365, 28, 28, 28, 90, 2…
$ compressive_strength <dbl> 79.99, 61.89, 40.27, 41.05, 44.30, 47.03, 43.70, …

Zmienna compressive_strength jest zmienną zależną, przewidywaną w tym zadaniu. W kilku przypadkach w tym zestawie danych, ta sama formuła betonu była testowana wielokrotnie. Wolimy nie uwzględniać tych replikowanych mieszanek jako pojedynczych punktów danych, ponieważ mogą one być rozmieszczone zarówno w zbiorze treningowym, jak i testowym. Może to sztucznie zawyżyć nasze szacunki wydajności.

Kod
concrete <- 
   concrete %>% 
   group_by(across(-compressive_strength)) %>% 
   summarize(compressive_strength = mean(compressive_strength),
             .groups = "drop")
nrow(concrete)
[1] 992

Podzielmy dane przy użyciu domyślnego stosunku 3:1 treningu do testu i ponownie wypróbujmy zestaw treningowy przy użyciu pięciu powtórzeń 10-krotnej walidacji krzyżowej:

Kod
set.seed(1501)
concrete_split <- initial_split(concrete, strata = compressive_strength)
concrete_train <- training(concrete_split)
concrete_test  <- testing(concrete_split)

set.seed(1502)
concrete_folds <- 
   vfold_cv(concrete_train, strata = compressive_strength, repeats = 5)

Niektóre modele (zwłaszcza sieci neuronowe, KNN i SVM) wymagają predyktorów, które zostały wyśrodkowane i przeskalowane, więc niektóre przepływy pracy modelu będą wymagały przepisów z tymi krokami przetwarzania wstępnego. Dla innych modeli, tradycyjne rozwinięcie modelu powierzchni odpowiedzi (tj. interakcje kwadratowe i dwukierunkowe) jest dobrym pomysłem. Dla tych celów tworzymy dwie receptury:

Kod
normalized_rec <- 
   recipe(compressive_strength ~ ., data = concrete_train) %>% 
   step_normalize(all_predictors()) 

poly_recipe <- 
   normalized_rec %>% 
   step_poly(all_predictors()) %>% 
   step_interact(~ all_predictors():all_predictors())

13.2 Określenie modeli

Teraz zdefiniujmy model, które chcemy przetestować:

Kod
library(rules)
library(baguette)

linear_reg_spec <- 
   linear_reg(penalty = tune(), mixture = tune()) %>% 
   set_engine("glmnet")

nnet_spec <- 
   mlp(hidden_units = tune(), penalty = tune(), epochs = tune()) %>% 
   set_engine("nnet", MaxNWts = 2600) %>% 
   set_mode("regression")

mars_spec <- 
   mars(prod_degree = tune()) %>%  #<- use GCV to choose terms
   set_engine("earth") %>% 
   set_mode("regression")

svm_r_spec <- 
   svm_rbf(cost = tune(), rbf_sigma = tune()) %>% 
   set_engine("kernlab") %>% 
   set_mode("regression")

svm_p_spec <- 
   svm_poly(cost = tune(), degree = tune()) %>% 
   set_engine("kernlab") %>% 
   set_mode("regression")

knn_spec <- 
   nearest_neighbor(neighbors = tune(), dist_power = tune(), weight_func = tune()) %>% 
   set_engine("kknn") %>% 
   set_mode("regression")

cart_spec <- 
   decision_tree(cost_complexity = tune(), min_n = tune()) %>% 
   set_engine("rpart") %>% 
   set_mode("regression")

bag_cart_spec <- 
   bag_tree() %>% 
   set_engine("rpart", times = 50L) %>% 
   set_mode("regression")

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

xgb_spec <- 
   boost_tree(tree_depth = tune(), learn_rate = tune(), loss_reduction = tune(), 
              min_n = tune(), sample_size = tune(), trees = tune()) %>% 
   set_engine("xgboost") %>% 
   set_mode("regression")

cubist_spec <- 
   cubist_rules(committees = tune(), neighbors = tune()) %>% 
   set_engine("Cubist") 

Autorzy w Khun i Johnson (2013) określa, że sieć neuronowa powinna mieć do 27 jednostek ukrytych w warstwie. Funkcja extract_parameter_set_dials() wyodrębnia zbiór parametrów, który modyfikujemy, aby miał prawidłowy zakres parametrów:

Kod
nnet_param <- 
   nnet_spec %>% 
   extract_parameter_set_dials() %>% 
   update(hidden_units = hidden_units(c(1, 27)))

W następnym kroku stworzymy zestawy przepływu pracy:

Kod
normalized <- 
   workflow_set(
      preproc = list(normalized = normalized_rec), 
      models = list(SVM_radial = svm_r_spec, SVM_poly = svm_p_spec, 
                    KNN = knn_spec, neural_network = nnet_spec)
   )
normalized
# A workflow set/tibble: 4 × 4
  wflow_id                  info             option    result    
  <chr>                     <list>           <list>    <list>    
1 normalized_SVM_radial     <tibble [1 × 4]> <opts[0]> <list [0]>
2 normalized_SVM_poly       <tibble [1 × 4]> <opts[0]> <list [0]>
3 normalized_KNN            <tibble [1 × 4]> <opts[0]> <list [0]>
4 normalized_neural_network <tibble [1 × 4]> <opts[0]> <list [0]>

Ponieważ zastosowaliśmy tylko jedna funkcję wstępnej obróbki danych (normalized_rec), to w podsumowaniu występują tylko kombinacje tego preprocesora i modeli. Kolumna wflow_id tworzona jest automatycznie, ale może być modyfikowana poprzez wywołanie mutate(). Kolumna info zawiera tibble z pewnymi identyfikatorami i obiektem przepływu pracy. Przepływ pracy może zostać wyodrębniony:

Kod
normalized %>% extract_workflow(id = "normalized_KNN")
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: nearest_neighbor()

── Preprocessor ────────────────────────────────────────────────────────────────
1 Recipe Step

• step_normalize()

── Model ───────────────────────────────────────────────────────────────────────
K-Nearest Neighbor Model Specification (regression)

Main Arguments:
  neighbors = tune()
  weight_func = tune()
  dist_power = tune()

Computational engine: kknn 

Kolumna option to miejsce na dowolne argumenty, których należy użyć, gdy oceniamy przepływ pracy. Na przykład, aby dodać obiekt parametrów sieci neuronowej:

Kod
normalized <- 
   normalized %>% 
   option_add(param_info = nnet_param, id = "normalized_neural_network")
normalized
# A workflow set/tibble: 4 × 4
  wflow_id                  info             option    result    
  <chr>                     <list>           <list>    <list>    
1 normalized_SVM_radial     <tibble [1 × 4]> <opts[0]> <list [0]>
2 normalized_SVM_poly       <tibble [1 × 4]> <opts[0]> <list [0]>
3 normalized_KNN            <tibble [1 × 4]> <opts[0]> <list [0]>
4 normalized_neural_network <tibble [1 × 4]> <opts[1]> <list [0]>

Kolumna result jest miejscem na wyjście funkcji dostrajania lub resamplingu. Dla innych modeli nieliniowych utwórzmy kolejny zestaw przepływów pracy:

Kod
model_vars <- 
   workflow_variables(outcomes = compressive_strength, 
                      predictors = everything())

no_pre_proc <- 
   workflow_set(
      preproc = list(simple = model_vars), 
      models = list(MARS = mars_spec, CART = cart_spec, CART_bagged = bag_cart_spec,
                    RF = rf_spec, boosting = xgb_spec, Cubist = cubist_spec)
   )
no_pre_proc
# A workflow set/tibble: 6 × 4
  wflow_id           info             option    result    
  <chr>              <list>           <list>    <list>    
1 simple_MARS        <tibble [1 × 4]> <opts[0]> <list [0]>
2 simple_CART        <tibble [1 × 4]> <opts[0]> <list [0]>
3 simple_CART_bagged <tibble [1 × 4]> <opts[0]> <list [0]>
4 simple_RF          <tibble [1 × 4]> <opts[0]> <list [0]>
5 simple_boosting    <tibble [1 × 4]> <opts[0]> <list [0]>
6 simple_Cubist      <tibble [1 × 4]> <opts[0]> <list [0]>

Na koniec składamy zestaw wykorzystujący warunki nieliniowe i interakcje z odpowiednimi modelami:

Kod
with_features <- 
   workflow_set(
      preproc = list(full_quad = poly_recipe), 
      models = list(linear_reg = linear_reg_spec, KNN = knn_spec)
   )

Te obiekty to tibble z dodatkową klasą workflow_set. Łączenie wierszy nie wpływa na stan zestawów, a wynik jest sam w sobie zestawem przepływów pracy:

Kod
all_workflows <- 
   bind_rows(no_pre_proc, normalized, with_features) %>% 
   # Make the workflow ID's a little more simple: 
   mutate(wflow_id = gsub("(simple_)|(normalized_)", "", wflow_id))
all_workflows
# A workflow set/tibble: 12 × 4
   wflow_id             info             option    result    
   <chr>                <list>           <list>    <list>    
 1 MARS                 <tibble [1 × 4]> <opts[0]> <list [0]>
 2 CART                 <tibble [1 × 4]> <opts[0]> <list [0]>
 3 CART_bagged          <tibble [1 × 4]> <opts[0]> <list [0]>
 4 RF                   <tibble [1 × 4]> <opts[0]> <list [0]>
 5 boosting             <tibble [1 × 4]> <opts[0]> <list [0]>
 6 Cubist               <tibble [1 × 4]> <opts[0]> <list [0]>
 7 SVM_radial           <tibble [1 × 4]> <opts[0]> <list [0]>
 8 SVM_poly             <tibble [1 × 4]> <opts[0]> <list [0]>
 9 KNN                  <tibble [1 × 4]> <opts[0]> <list [0]>
10 neural_network       <tibble [1 × 4]> <opts[1]> <list [0]>
11 full_quad_linear_reg <tibble [1 × 4]> <opts[0]> <list [0]>
12 full_quad_KNN        <tibble [1 × 4]> <opts[0]> <list [0]>

13.3 Tuning modeli

Prawie wszystkie modele ujęte w all_workflows zawierają parametry dostrajania. Aby ocenić ich wydajność, możemy użyć standardowych funkcji strojenia lub resamplingu (np. tune_grid() i tak dalej). Funkcja workflow_map() zastosuje tę samą funkcję do wszystkich przepływów w zestawie; domyślnie jest to tune_grid().

Dla tego przykładu, wyszukiwanie w oparciu o siatkę jest stosowane do każdego przepływu pracy, stosując jednocześnie 25 różnych kandydatów na parametry. Istnieje zestaw wspólnych opcji do wykorzystania przy każdym wykonaniu tune_grid(). Na przykład, w poniższym kodzie użyjemy tego samego próbkowania i obiektów kontrolnych dla każdego przepływu pracy, wraz z rozmiarem siatki równym 25. Funkcja workflow_map() posiada dodatkowy argument o nazwie seed, który służy do zapewnienia, że każde wykonanie tune_grid() zużywa tych samych liczb losowych.

Kod
library(doParallel)
cl <- makeCluster(4)
registerDoParallel(cl)

grid_ctrl <-
   control_grid(
      save_pred = TRUE,
      parallel_over = "everything",
      save_workflow = TRUE
   )

grid_results <-
   all_workflows %>%
   workflow_map(
      seed = 1503,
      resamples = concrete_folds,
      grid = 25,
      control = grid_ctrl
   )

W podsumowaniu widać, że kolumny option i result zostały zaktualizowane:

# A workflow set/tibble: 12 × 4
   wflow_id             info             option    result   
   <chr>                <list>           <list>    <list>   
 1 MARS                 <tibble [1 × 4]> <opts[3]> <tune[+]>
 2 CART                 <tibble [1 × 4]> <opts[3]> <tune[+]>
 3 CART_bagged          <tibble [1 × 4]> <opts[3]> <rsmp[+]>
 4 RF                   <tibble [1 × 4]> <opts[3]> <tune[+]>
 5 boosting             <tibble [1 × 4]> <opts[3]> <tune[+]>
 6 Cubist               <tibble [1 × 4]> <opts[3]> <tune[+]>
 7 SVM_radial           <tibble [1 × 4]> <opts[3]> <tune[+]>
 8 SVM_poly             <tibble [1 × 4]> <opts[3]> <tune[+]>
 9 KNN                  <tibble [1 × 4]> <opts[3]> <tune[+]>
10 neural_network       <tibble [1 × 4]> <opts[4]> <tune[+]>
11 full_quad_linear_reg <tibble [1 × 4]> <opts[3]> <tune[+]>
12 full_quad_KNN        <tibble [1 × 4]> <opts[3]> <tune[+]>

Kolumna option zawiera teraz wszystkie opcje, których użyliśmy w wywołaniu workflow_map(). W kolumnach result, notacje tune[+] i rsmp[+] oznaczają, że obiekt nie miał żadnych problemów w procesie optymalizacji. Wartość taka jak tune[x] pojawia się, gdy wszystkie modele z jakiegoś powodu zawiodły.

Istnieje kilka wygodnych funkcji do badania wyników, takich jak grid_results. Funkcja rank_results() uporządkuje modele według wybranej metryki wydajności. Domyślnie używa ona pierwszej metryki w zestawie metryk (w tym przypadku RMSE). Przefiltrujmy wyniki, aby analizować tylko na RMSE:

Kod
grid_results %>% 
   rank_results() %>% 
   filter(.metric == "rmse") %>% 
   select(model, .config, rmse = mean, rank)
# A tibble: 252 × 4
   model        .config                rmse  rank
   <chr>        <chr>                 <dbl> <int>
 1 boost_tree   Preprocessor1_Model04  4.25     1
 2 boost_tree   Preprocessor1_Model06  4.29     2
 3 boost_tree   Preprocessor1_Model13  4.31     3
 4 boost_tree   Preprocessor1_Model14  4.39     4
 5 boost_tree   Preprocessor1_Model16  4.46     5
 6 boost_tree   Preprocessor1_Model03  4.47     6
 7 boost_tree   Preprocessor1_Model15  4.48     7
 8 boost_tree   Preprocessor1_Model05  4.55     8
 9 boost_tree   Preprocessor1_Model20  4.71     9
10 cubist_rules Preprocessor1_Model24  4.71    10
# ℹ 242 more rows

Domyślnie funkcja szereguje wszystkie zestawy kandydatów, dlatego ten sam model może pojawić się wielokrotnie na wyjściu. Opcja select_best może być użyta do uszeregowania modeli przy użyciu najlepszej kombinacji dostrajania parametrów. Metoda autoplot() tworzy wykresy rankingowy; posiada ona również argument select_best. Wykres na Rys. 13.1 wizualizuje najlepsze wyniki dla każdego modelu i jest generowany za pomocą:

Kod
autoplot(
   grid_results,
   rank_metric = "rmse",  # <- how to order models
   metric = "rmse",       # <- which metric to visualize
   select_best = TRUE     # <- one point per workflow
) +
   geom_text(aes(y = mean - 1/2, label = wflow_id), angle = 90, hjust = 1) +
   lims(y = c(3, 9.5)) +
   theme(legend.position = "none")
Rys. 13.1: Oszacowany RMSE (i przybliżone przedziały ufności) dla najlepszej konfiguracji modelu w każdym przepływie pracy

W przypadku, gdy chcesz zobaczyć wyniki dostrajania parametrów dla konkretnego modelu, tak jak na Rys. 13.2, argument id może przyjąć pojedynczą wartość z kolumny wflow_id dla którego modelu ma być wykreślony:

Kod
autoplot(grid_results, id = "Cubist", metric = "rmse")
Rys. 13.2: Wizualizacja RMSE w kontekście konfiguracji hiperparametrów modelu Cubist

Funkcje collect_metrics() i collect_predictions() również pozwalają przejrzenie wyników optymalizacji.

W powyższym procesie optymalizacji przeuczono 12600 modeli, co zajęło około 2 godzin przy wykorzystaniu 4 rdzeni procesora. Pokazuje to, że zagadnienie tuningu nawet kilku kandydackich modeli zajmuje sporo czasu.

13.4 Efektywna filtracja modeli

Jedną z metod efektywnego przesiewania dużego zbioru modeli jest zastosowanie podejścia wyścigowego opisanego wcześniej. Mając zestaw przepływów pracy, możemy użyć funkcji workflow_map() do podejścia wyścigowego.

Kod
race_ctrl <-
   control_race(
      save_pred = TRUE,
      parallel_over = "everything",
      save_workflow = TRUE
   )

race_results <-
   all_workflows %>%
   workflow_map(
      "tune_race_anova",
      seed = 1503,
      resamples = concrete_folds,
      grid = 25,
      control = race_ctrl
   )

Nowy obiekt wygląda bardzo podobnie, choć elementy kolumny wyników wykazują wartość race[+], co wskazuje na inny typ obiektu:

Kod
race_results <- readRDS("models/race_results.rds")
race_results
# A workflow set/tibble: 12 × 4
   wflow_id             info             option    result   
   <chr>                <list>           <list>    <list>   
 1 MARS                 <tibble [1 × 4]> <opts[3]> <race[+]>
 2 CART                 <tibble [1 × 4]> <opts[3]> <race[+]>
 3 CART_bagged          <tibble [1 × 4]> <opts[3]> <rsmp[+]>
 4 RF                   <tibble [1 × 4]> <opts[3]> <race[+]>
 5 boosting             <tibble [1 × 4]> <opts[3]> <race[+]>
 6 Cubist               <tibble [1 × 4]> <opts[3]> <race[+]>
 7 SVM_radial           <tibble [1 × 4]> <opts[3]> <race[+]>
 8 SVM_poly             <tibble [1 × 4]> <opts[3]> <race[+]>
 9 KNN                  <tibble [1 × 4]> <opts[3]> <race[+]>
10 neural_network       <tibble [1 × 4]> <opts[4]> <race[+]>
11 full_quad_linear_reg <tibble [1 × 4]> <opts[3]> <race[+]>
12 full_quad_KNN        <tibble [1 × 4]> <opts[3]> <race[+]>
Kod
autoplot(
   race_results,
   rank_metric = "rmse",  
   metric = "rmse",       
   select_best = TRUE    
) +
   geom_text(aes(y = mean - 1/2, label = wflow_id), angle = 90, hjust = 1) +
   lims(y = c(3.0, 9.5)) +
   theme(legend.position = "none")
Rys. 13.3: Oszacowane RMSE (i przybliżone przedziały ufności) dla najlepszej konfiguracji modelu w poszukiwaniu za pomocą metody wyścigowej.

Podejście wyścigowe oszacowało łącznie 1050 modeli, 8,33% z pełnego zestawu 12600 modeli w pełnej siatce. W rezultacie podejście wyścigowe trwało nieco ponad 17 min., więc było 7-krotnie szybsze1.

1 Wartości te będą zależały od sprzętu na jakim wykonuje się obliczenia

Na ile zbliżone wyniki otrzymaliśmy stosując obie metody tuningu?

Kod
matched_results <- 
   rank_results(race_results, select_best = TRUE) %>% 
   select(wflow_id, .metric, race = mean, config_race = .config) %>% 
   inner_join(
      rank_results(grid_results, select_best = TRUE) %>% 
         select(wflow_id, .metric, complete = mean, 
                config_complete = .config, model),
      by = c("wflow_id", ".metric"),
   ) %>%  
   filter(.metric == "rmse")

library(ggrepel)

matched_results %>% 
   ggplot(aes(x = complete, y = race)) + 
   geom_abline(lty = 3) + 
   geom_point() + 
   geom_text_repel(aes(label = model)) +
   coord_obs_pred() + 
   labs(x = "Complete Grid RMSE", y = "Racing RMSE") 

Podczas gdy podejście wyścigowe wybrało te same parametry kandydata co kompletna siatka tylko dla 41,67% modeli, metryki wydajności modeli wybranych przez wyścig były prawie równe. Korelacja wartości RMSE wyniosła 0,968, a korelacja rangowa 0,951. Wskazuje to, że w obrębie modelu istniało wiele kombinacji dostrajania parametrów, które dawały niemal identyczne wyniki.

13.5 Finalizacja modelu

Podobnie do tego, co pokazaliśmy w poprzednich rozdziałach, proces wyboru ostatecznego modelu i dopasowania go na zbiorze treningowym jest prosty. Pierwszym krokiem jest wybranie zbioru treningowego do sfinalizowania. Ponieważ model boosted tree działał dobrze, wyodrębnimy go ze zbioru, zaktualizujemy parametry o numerycznie najlepsze ustawienia i dopasujemy do zbioru treningowego. W przypadku gdy mamy wątpliwości dotyczące siatki hiperparametrów dobranych podczas filtrowania modeli, np. że pomija ona ważne kombinacje, możemy zastosować do wybranego modelu metody finetune przedstawione w poprzednim rozdziale.

Kod
best_results <- 
   race_results %>% 
   extract_workflow_set_result("boosting") %>% 
   select_best(metric = "rmse")
best_results
# A tibble: 1 × 7
  trees min_n tree_depth learn_rate loss_reduction sample_size .config          
  <int> <int>      <int>      <dbl>          <dbl>       <dbl> <chr>            
1  1957     8          7     0.0756    0.000000145       0.679 Preprocessor1_Mo…
Kod
boosting_test_results <- 
   race_results %>% 
   extract_workflow("boosting") %>% 
   finalize_workflow(best_results) %>% 
   last_fit(split = concrete_split)

Wyniki metryki dla zbioru testowego oraz wizualizację predykcji możemy zobaczyć stosując.

Kod
collect_metrics(boosting_test_results)
# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       3.41  Preprocessor1_Model1
2 rsq     standard       0.954 Preprocessor1_Model1
Kod
boosting_test_results %>% 
   collect_predictions() %>% 
   ggplot(aes(x = compressive_strength, y = .pred)) + 
   geom_abline(color = "gray50", lty = 2) + 
   geom_point(alpha = 0.5) + 
   coord_obs_pred() + 
   labs(x = "observed", y = "predicted")
Rys. 13.4: Porównanie wartości obserwowanych i przewidywanych z modelu