Close
Enter your
text here

Machine Learning

Ako zrýchliť XGBoost v R: Benchmark trénovania modelu na CPU vs. GPU 150 150 cleandata

Ako zrýchliť XGBoost v R: Benchmark trénovania modelu na CPU vs. GPU
cpu-vs-gpu-blog

Ako rýchlo dokážete vytrénovať model na veľkých datasetoch? Toto je otázka, ktorú si kladie mnoho analytikov a data scientistov, keď začínajú pracovať s pokročilými algoritmami strojového učenia, ako je XGBoost. V tomto článku si ukážeme, ako vytvoriť klasifikačný model pomocou XGBoost v R a porovnáme si dĺžku výpočtu pri použití CPU a GPU.

Štandardná verzia knižnice XGBoost na CRAN však nepodporuje parameter device = “cuda”. Preto je potrebné nainštalovať špeciálnu verziu XGBoostu – buď skompilovať z repozitára na GitHube, alebo využiť poskytovaný inštalačný balík (napríklad pre Windows takto):

Inštalácia xgboost 
options(timeout = 1200)
xgboost_url <- "https://s3-us-west-2.amazonaws.com/xgboost-nightly-builds/release_2.0.0/xgboost_r_gpu_win64_82d846bbeb83c652a0b1dff0e3519e67569c4a3d.tar.gz"
install.packages(xgboost_url, repos = NULL, type = "source")

CPU vs GPU

Na veľkých datasetoch môže trénovanie modelov trvať hodiny, ak nie dni, najmä na bežných procesoroch (CPU). Moderné grafické karty (GPU) od spoločnosti Nvidia však umožňujú paralelné spracovanie dát, čo dramaticky zrýchľuje niektoré výpočtové operácie.

CPU je ako pracovník, ktorý rieši úlohy sekvenčne, zatiaľ čo GPU je tím tisícov pracovníkov, ktorí dokážu riešiť mnohé úlohy naraz.

Ešte predtým, než prejdeme k modelovaniu, môžeme si zistiť, aký hardvér vlastne používame.

AMD Ryzen 9 7900X je vysoko výkonný desktopový procesor s 12 fyzickými jadrami (24 vláknami), ktorý je všeobecne považovaný za high-end/výkonný segment pre multivláknové operácie (aktuálne je dostupná už novšia generácia AMD procesorov). K nemu je 67.8 GB RAM, čo dáva dostatočnú rezervu na spracovanie stredne veľkých až väčších datasetov. Grafická karta Nvidia RTX 4060Ti s 8 GB VRAM sa radí do strednej až vyššej strednej triedy, no pri veľmi rozsiahlych datasetoch môže byť 8 GB VRAM limitujúce. Napriek tomu dokáže pri bežných projektoch a primerane veľkých dátach výrazne zrýchliť trénovanie XGBoost modelov v porovnaní s CPU. Viac info nájdete v staršom blogu.

1. Načítanie knižníc a príprava prostredia

Na začiatku si načítame potrebné balíky. Je dobré nastaviť aj niektoré globálne možnosti, napríklad options(scipen = 999), aby sme obmedzili vedecký zápis čísel.

Knižnice 
options(scipen = 999)
library(data.table)
library(readr)
library(xgboost)
library(tidymodels)
library(dplyr)
library(stringr)
library(doParallel)
library(tictoc)

data.table nám umožní rýchle nahratie .csv. tidymodels je framework pre tvorbu a ladenie modelov. xgboost je samotná implementácia gradient boosting algoritmov. tictoc nám zmeria čas operácie (funkcie tic() a toc())

2. Načítanie, čistenie a rozdelenie dát

Začneme s dvoma datasetmi: train.csv a test.csv. Oba súbory sú dostupné na Kaggle.

Dátový stĺpec id odstránime (neobsahuje relevantné informácie pre model). Premennú class transformujeme tak, aby hodnota “e” mala label 0, a “p” label 1. Prázdne stringy (textové premenné) nahradíme za “missing”, pretože chceme explicitne ošetriť chýbajúce hodnoty.

Príprava dát 
## data prep ----
train <- data.table::fread("data/train.csv") |> 
  select(-id) |> 
  mutate(class = factor(class, levels = c("p", "e"), labels = c(1, 0))) |> 
  mutate(across(where(is.character), ~ if_else(. == "", "missing", .)))

dim(train)
[1] 3116945      21

Finálny trénovací dataset má 21 stĺpcov vrátane triedy, ktorú budeme predikovať a 3 116 945 záznamov.

Nasleduje klasické rozdelenie pomocou initial_split z rsample: prop = 0.7 znamená, že 70 % záznamov ide do trénovacej množiny a 30 % do testovacej, strata = class zachováva približný pomer tried v oboch množinách.

Pri ladení modelu použijeme 10-násobnú cross validáciu.

Rozdelenie dát 
set.seed(123)
split <- initial_split(train, prop = 0.7, strata = class)

train_set <- training(split)
test_set <- testing(split)

set.seed(123)
folds <- vfold_cv(train_set, strata = class)

3. Pre-Processing

V recipe() si definujeme jednotlivé kroky v tzv. recipe pipeline.

step_other() s threshold = 0.005 zlúči veľmi vzácne kategórie (s podielom < 0.5 %) do kategórie “other”. step_string2factor() prevedie stringové atribúty na faktory. step_dummy() vytvorí dummy premenné z faktorov.

Nastavenie pre-processingu 
recipe <- recipe(class ~ ., data = train_set) |> 
  step_other(all_nominal_predictors(), threshold = 0.005) |> 
  step_string2factor(all_nominal_predictors()) |>
  step_dummy(all_nominal_predictors())

4. Definícia XGBoost modelov (CPU a GPU)

4.1. Základné XGBoost parametre

  • trees = 1000: Počet stromov.
  • tree_depth = tune(): Hĺbka stromu, ktorú budeme ladiť.
  • min_n = tune(): Minimálny počet záznamov v koncovom uzle stromu (v XGBoost často min_child_weight).
  • loss_reduction = tune(): (v XGBoost označované ako gamma) – ovplyvňuje, o koľko sa musí znížiť chyba, aby došlo k deleniu uzla.
  • sample_size = tune(): (v XGBoost subsample) – podiel vzoriek použitých pre každý strom.
  • mtry = tune(): (v XGBoost colsample_bynode) – podiel stĺpcov vyberaných pri každom delení.
  • learn_rate = tune(): (v XGBoost eta) – reguluje rýchlosť učenia.

Rozdiel medzi CPU a GPU modelmi: Zásadné je device = “cuda” a tree_method = “hist”. Týmito parametrami vravíme xgboostu, že chceme použiť GPU. Nastavujeme nthread = 12, aby sme využili všetky dostupné jadrá CPU.

4.2. Model pre CPU

Definovanie CPU modelu 
xgb_cpu <- boost_tree(
  trees = 1000,
  tree_depth = tune(), min_n = tune(),
  loss_reduction = tune(),                     ## first three: model complexity
  sample_size = tune(), mtry = tune(),         ## randomness
  learn_rate = tune()                          ## step size
) |> 
  set_engine("xgboost", nthread = 12) |> 
  set_mode("classification")

4.3. Model pre GPU

Definovanie GPU modelu 
xgb_gpu <- boost_tree(
  trees = 1000,
  tree_depth = tune(), min_n = tune(),
  loss_reduction = tune(),                     ## first three: model complexity
  sample_size = tune(), mtry = tune(),         ## randomness
  learn_rate = tune()                          ## step size
) |> 
  set_engine("xgboost", device = "cuda", tree_method = "hist", nthread = 12) |> 
  set_mode("classification")

5. Workflows

V tidymodels je workflow elegantný spôsob, ako zlúčiť recept a model:

Vytvorenie workflows 
xgb_cpu_wf <- workflow() |> 
  add_recipe(recipe) |> 
  add_model(xgb_cpu)

xgb_gpu_wf <- workflow() |> 
  add_recipe(recipe) |> 
  add_model(xgb_gpu)

Aby sme mohli ladiť hyperparametre, vytvoríme si ešte bayes_param objekt. Použijeme finalize(mtry(), train_set), aby sme odhadli maximálnu hodnotu mtry na základe počtu stĺpcov.

Extrakcia hyperparametrov 
bayes_param <- xgb_gpu_wf %>% 
  extract_parameter_set_dials() %>% 
  update(mtry = finalize(mtry(), train_set))

6. Bayesovské ladenie a meranie času

Ladenie bude prebiehať v 6 iteráciách (parametrizovaných v iter = 6). Ide len o ukážku, v reálnom projekte by sme spravidla volili viac iterácií (10, 20, či viac), aby sme dôkladnejšie preskúmali priestor parametrov.

6.1. CPU – tune_bayes

Tuning CPU modelu 
tic()
set.seed(234)
cpu_xgb_res <- tune_bayes(
  xgb_cpu_wf,
  param_info = bayes_param,
  resamples = folds, 
  iter = 6
)
toc()
# 24278.11 sec elapsed

6.2. GPU – tune_bayes

Tuning GPU modelu 
tic()
set.seed(345)
gpu_xgb_res <- tune_bayes(
  xgb_gpu_wf,
  param_info = bayes_param,
  resamples = folds, 
  iter = 6
)
toc()
# 7345.68 sec elapsed

Porovnanie času výpočtu

  • CPU: ~ 6 hodín 45 minút
  • GPU: ~ 2 hodiny

GPU model je teda približne 3x rýchlejší, čo predstavuje významnú úsporu času pri práci na veľkých datasetoch. V praxi ešte záleží na type GPU, verzii CUDA, atď.

7. Výber najlepších parametrov a finálny model

Po skončení hyperparameter tuningu si môžeme pozrieť metriky cez collect_metrics(), a vybrať najlepšie nastavenia napr. na základe metriky roc_auc.

Výsledné CPU metriky 
metrics_cpu <- collect_metrics(cpu_xgb_res)
best_cpu <- select_best(cpu_xgb_res, metric = "roc_auc")

cpu_roc_auc <- metrics_cpu |> 
  filter(.metric == "roc_auc") |> 
  slice_max(mean) |> 
  pull(mean)

Pre CPU model je výsledná hodnota 0.9971878.

Výsledné GPU metriky 
metrics_gpu <- collect_metrics(gpu_xgb_res)
best_gpu <- select_best(gpu_xgb_res, metric = "roc_auc")

gpu_roc_auc <- metrics_gpu |> 
  filter(.metric == "roc_auc") |> 
  slice_max(mean) |> 
  pull(mean)

Pre GPU je to 0.9972345.

GPU model nielen že trénuje rýchlejšie, ale dosahuje aj lepšiu presnosť. Je to však veľmi individuálne, závisí od datasetu i parametrov.

Zhrnutie

Čo sme videli?

  • XGBoost sa dá používať na CPU aj GPU pomocou jednoduchých parametrov (device = “cuda”, tree_method = “hist”).
  • Čas výpočtu môže klesnúť až na zlomok, ak máme vhodnú GPU a väčší dataset.
  • Pri ladení parametrov sme použili Bayesovskú optimalizáciu, ktorá je síce efektívnejšia vo vyhľadávaní správnych hyperparametrov, ale je aj výpočtovo náročnejšia než jednoduchšie metódy.

Praktické tipy a nuansy

  • Pamäť GPU môže byť limitujúca. Veľké datasety nemusia vojsť do GPU pamäte.
  • Overheads: Pri relatívne malých datasetoch môže byť efekt GPU nižší alebo dokonca negatívny pre “overheads” spojené s kopírovaním dát medzi CPU a GPU.
  • Data pre-processing: Prechod na GPU je nápomocný hlavne pri trénovaní modelu, ale pre-processing zostáva často na CPU. Ak máte rozsiahlu prípravu dát, toto stále môže byť “bottleneck”. Aj z tohto dôvodu ponechávame nthreads na dostatočne vysokej hodnote.
  • Viac iterácií v Bayesovskej optimalizácii obvykle pomôže presnejšie nájsť najlepšie hyperparametre.

Záver

Modelovanie s XGBoost algoritmom na GPU môže drasticky skrátiť čas trénovania oproti CPU, najmä pri veľkých datasetoch. Ak je to aj Váš prípad a navyše pravidelne pracujete so zložitými modelmi, GPU je dlhodobou investíciou, ktorá šetrí čas a zvyšuje efektivitu práce. Pre menšie projekty však bude postačujúci aj výkon moderného CPU.

Forecasting so strojovým učením 150 150 cleandata

Forecasting so strojovým učením

V tomto blogu sa zameriavam na moje prvé skúsenosti s účasťou v forecast súťaži, kde som sa rozhodol otestovať nové metódy. Namiesto tradičných štatistických prístupov som sa sústredil na použitie modelov strojového učenia (ML). Hlavným cieľom bolo preskúmať, ako môžu moderné ML metódy zlepšiť presnosť predikcií v porovnaní s bežne používanými štatistickými technikami.

kaggle_forecast_zhrnutie

Moje prvé skúsenosti a lekcie z forecast súťaže od rohlik.cz

Nedávno som sa po prvýkrát zúčastnil forecastovacej súťaže na Kaggle. Organizovala ju spoločnosť Rohlik.cz. Aj keď môj výsledok bol pre mňa osobne sklamaním, keďže som skončil v 27. percentile, súťaž mi priniesla množstvo cenných skúseností a ponaučení, ktoré stoja za zmienku.

Prvé skúšanie ML na forecasting

Táto súťaž bola mojou prvou príležitosťou vyskúšať forecasting pomocou strojového učenia. Doteraz som využíval najmä tradičné štatistické metódy ako ARIMA a ETS. Tentokrát som sa rozhodol pre strojové učenie z nasledujúcich dôvodov:

  • Cieľom bolo forecastovať denné objednávky. S tradičnými metódami mám dobrú skúsenosť, pokiaľ sa jedná napr. o mesačnú agregáciu a mám dostupné historické dáta za dostatočne dlhé obdobie (niekoľko rokov). V tomto prípade sme mali dostupné údaje za relatívne krátke historické obdobie.
  • Strojové učenie bolo víťaznou stratégiou v mnohých predchádzajúcich súťažiach ako napr. M4, M5 alebo v súťaži organizovanej spoločnosťou Intermarché. Preto považujem za dôležité mať praktickú skúsenosť aj s týmto spôsobom forecastovania.
  • Už dlhšie som sa chystal vyskúšať knižnicu modeltime a ďalšie, k nej doplnkové knižnice.
  • Praktická skúsenosť s novou metódou mi pomôže aj v pracovnej sfére.

Priebeh a výsledok

Na súťaži som pracoval vo voľnom čase, čo znamenalo, že som musel prioritizovať, čomu sa chcem venovať.

Základ by ideálne bol rozsiahly a kvalitný feature engineering. Spätne musím povedať, že som mu mal venovať viac času. Rozdiely v konečnom hodnotení boli relatívne malé a verím, že práve lepšie features by ma vedeli posunúť na vyššiu priečku. Avšak hlavným cieľom bolo vyskúšať nové knižnice a metódy, takže som venoval až neprimerane veľa času (ak by sa jednalo o reálny projekt) ich implementácii.

Hodnotená metrika bola MAPE. Nie je to ideálna metrika, keďže jednotlivé sklady, pre ktoré sa robili forecasty, mali veľmi rozdielny objem objednávok (od dolných tisícov po vyse 11 tisíc) a teda 10% odchýlka v jednom sklade mohla byť 200 objednávok a 4% odchýlka v inom 400 objednávok. Moja výsledná hodnota MAPE bola 0,0489, pričom v priebežnom public leaderboarde bola 0.0396 a v cross-validácii 0,0331. Napriek tomuto zhoršeniu sa mi priečka posunula len o 8 miest dole.

Problém so súťažami na Kaggle je, že mnoho súťažiacich skopíruje dobré uverejnené riešenia iných a nahrajú výsledky. Na grafe nižšie je efekt pekne viditeľný v intervale 4,6% až 4,9% a bol často kritizovaný v diskusiách. Nie je nič zlé na inšpirovaní sa u iných, avšak skopírovanie celého kódu s kozmetickými úpravami je neetické. V grafe som zvýraznil interval, v ktorom sa nachádzalo moje riešenie. Objektívne však menej ako 5% odchýlku na denných predajoch pre 7 skladov na nasledujúcich 60 dní považujem za dobrý výsledok. Graf tiež dobre ilustruje ako blízko boli riešenia pri sebe. Šírka intervalov je 0,25 percenta.

Distribúcia výsledných MAPE v private leaderboarde

Testovanie rôznych prístupov

Zvolil som dva rôzne prístupy, ktoré som chcel otestovať.

Prvý bol tradičnejší workflow – feature engineering, zadefinovanie modelu, tuning hyperparametrov s cross-validáciou a následný testing na testovacej sade dát. Použil som tri rôzne algoritmy: XGBoost, LightGBM a Random Forest. Každý z týchto modelov má svoje vlastné silné stránky, no výsledky ukázali, že v priemere najlepší bol LightGBM, mierne horší XGBoost a po nich s väčším prepadom v presnosti Random Forest. Keďže sa však výsledky jednotlivých modelov mierne líšili podľa skladu, pre ktorý bol forecast robený, rozhodol som sa použiť ensemble prístup s meta learner modelom (použitý algoritmus bol GLMnet). Tento model použil výsledky z troch algoritmov plus dve premenné – sklad a deň v týždni. Je celkom bežné, že aj obyčajné spriemerovanie nezávislých forecastov zlepší výslednú presnosť.

Druhý prístup bol rekurzívny ML model. Ten obsahuje všetky kroky z predchádzajúceho workflowu, ale pridáva sa funkcia, ktorá dopočítava niektoré hodnoty z predchádzajúcich dní. Napr. predaje v predchádzajúci deň, 7 dní dozadu, dva týždne dozadu atď. Tento proces je pomalší, lebo vyžaduje, aby sa robil forecast po jednom dni, ten sa následne “doplní” do dátového setu a vstupuje do nasledujúceho forecastu už ako prediktor. Jedná sa teda o iteratívny proces. Keďže tento model postupne predikoval hodnoty pre jednotlivé časové obdobia na základe predchádzajúcich predikcií, jeho presnosť závisí od presnosti predchádzajúcich predikcií, čo prináša určité riziká.

Jedným špecifikom forecastov časových radov je, že cross-validácia (CV) by mala rešpektovať kontinuitu času. Time series cross-validation môže mať niekoľko podôb, napr. sliding alebo expanding window, testovacie dáta môžu priamo nadväzovať na trénovacie alebo môžeme aplikovať lag. Spoločná myšlienka za nimi však je, že pri trénovaní modelov pre predpovedanie časových radov by v prípade klasických metód ako k-fold CV dochádzalo k tzv. data leakage (trénovanie modelu na dátach z budúcnosti a teda na informáciách, ktoré by model nemal mať). Vačšina riešení, ktoré sú zverejnené sa však tejto téme nevenuje a rovno uvádza hodnoty hyperpametrov, ktoré boli aplikované.

Dôležitosť feature engineeringu

Jedným z najzaujímavejších zistení bolo, že viacero z top riešení v súťaži používalo len jeden model – LightGBM. Najlepšie riešenie nemá príliš dobrú dokumentáciu, avšak toto riešenie skončilo na 2. mieste a potvrdilo, že niekedy menej je viac. Tento úspech ukazuje, že najdôležitejším krokom je správny feature engineering. V niektorých prípadoch je lepšie investovať viac času do správnej prípravy dát, než do komplikovaných modelov, ktoré sú v podstate „blackboxom“. V tejto súťaži sa navyše nedalo príliš spoliehať na priebežné výsledky v public leaderboarde, keďže validačný set dát mal len 31 % z už tak krátkeho forecastovaného obdobia. Namiesto toho bolo lepšie dôverovať výsledkom vlastnej cross-validácie.

Záver

Aj keď som nedosiahol výsledky v aké som dúfal (byť aspoň v top 20 percent), súťaž mi dala možnosť vyskúšať si nové techniky, zdokonaliť svoje schopnosti v oblasti strojového učenia a získať cenné skúsenosti. Ak sa zúčastňujete podobných súťaží, nezabudnite, že aj keď váš model nedosiahne top skóre, stále sa môžete veľa naučiť a získať inšpiráciu z vyššie umiestnených riešení.

Machine learning v jazyku R – Odhad cien bytov 150 150 cleandata

Machine learning v jazyku R – Odhad cien bytov

Machine Learning (Strojové učenie) je podmnožina umelej inteligencie, ktorá sa zameriava na vývoj algoritmov schopných učiť sa a robiť predikcie na základe dát. V jazyku R existuje množstvo nástrojov a frameworkov, ktoré uľahčujú vytváranie a tréning ML modelov. Medzi najpopulárnejšie patrí tidymodels, ktorý poskytujú robustné a efektívne prostriedky na tvorbu rôznych modelov a predikciu. V tomto článku sa zameriam na vytvorenie modelu na predikovanie cien bytov, pričom priblížim aj tému explainable ML na lepšie porozumenie výsledkom modelu.

Vytvorenie ML Modelu v jazyku R

Úvod

V tomto blogovom príspevku prejdem procesom prípravy a tvorby ML modelu, ktorý bude na základe vložených parametrov predikovať cenu bytu. Dáta boli vopred nachystané, ich čistenie a manipuláciu si môžete pozrieť v tomto článku.

Explainable ML

V rámci strojového učenia predstavuje koncept “explainable/interpretable ML” (vysvetliteľné strojové učenie) kľúčový posun smerom k transparentnosti a interpretovateľnosti modelov. Napriek tomu, že modely ako XGBoost ponúkajú výnimočnú prediktívnu silu, ich interné rozhodovacie procesy môžu byť zložité a nejasné, a teda môžu byť tiež vnímané ako “black box”, čo znamená, že nie je ľahké pochopiť, ako sa k predikciám dospelo. Toto vnímanie komplikuje dôveru a akceptáciu modelov v kritických aplikáciách, kde je potrebné pochopenie dôvodov za predikciami.

Vysvetliteľné strojové učenie sa snaží preklenúť túto priepasť poskytujúc nástroje a metódy na objasnenie, ako modely dospeli k svojim rozhodnutiam. Knižnice ako “vip” v jazyku R a metódy ako SHAP (Shapley Additive exPlanations) hodnoty umožňujú analytikom a vývojárom lepšie pochopiť príspevky jednotlivých funkcií k výslednému predikovanému výstupu. Táto schopnosť detailne rozložiť predikčný proces umožňuje nielen hlbšiu analýzu a optimalizáciu modelov, ale tiež zvyšuje transparentnosť a dôveru zo strany koncových užívateľov.

Vysvetliteľné strojové učenie. Zdroj: Lundberg, Scott & Erion, Gabriel & Chen, Hugh & DeGrave, Alex & Prutkin, Jordan & Nair, Bala & Katz, Ronit & Himmelfarb, Jonathan & Bansal, Nisha & Lee, Su-In. (2019). Explainable AI for Trees: From Local Explanations to Global Understanding.

Dôležitosť explainable/interpretable ML naberá na váhe najmä v sektoroch, kde sú dôsledky rozhodnutí založených na predikciách modelu vysoké, ako sú zdravotníctvo, financie alebo právo. V týchto oblastiach je kľúčové, aby boli modely nielen presné, ale aj ich rozhodnutia pochopiteľné a spravodlivé. Vysvetliteľné strojové učenie tak stojí v centre úsilia o vytvorenie technológií, ktoré sú nielen inteligentné, ale aj zrozumiteľné a etické.

Výborná kniha na túto tému je napríklad Interpretable Machine Learning: A Guide For Making Black Box Models Explainable, ktorá je veľmi zrorumiteľne napísaná a obsahuje aj odkazy na knižnice v R aj Pythone. Navyše je dostupná aj online. Náročnejšia publikácia je potom Explainable AI for Trees: From Local Explanations to Global Understanding, ktorá ide viac do detailu a “matematiky” na pozadí.

Načítanie knižníc a dát

Klasicky začíname s načítaním knižníc. Aj tento krát použijeme tidyverse – súbor knižníc navrhnutých na prácu s dátami v R, ktorý zahŕňa napr. ggplot2 pre vizualizáciu, dplyr pre manipuláciu s dátami, tidyr pre úpravu tvaru dát, a iné. Opäť sa objavuje knižnica sf, keďže dáta obsahujú priestorové objekty. Novikou v tomto blogu je tidymodels. Je to framework pre modelovanie a strojové učenie, ktorý poskytuje koherentnú súpravu knižníc na predspracovanie dát, rozdelenie dát, cross-validáciu, výber modelu atď. Je navrhnutý tak, aby bol v súlade s princípmi tidyverse a umožňoval ľahkú integráciu s inými nástrojmi z tohto ekosystému. Dve nasledujúce knižnice pomáhajú s konceptom “explainable/interpretable ML”. Knižnica vip je určená na vizualizáciu dôležitosti premenných v rôznych modeloch strojového učenia. Umožňuje ľahko identifikovať, ktoré premenné majú najväčší vplyv na predikcie modelu, čo je kľúčové pre interpretáciu modelu a pochopenie dát. Knižnica shapviz je špecificky navrhnutá na výpočet a vizualizáciu SHAP hodnôt pre modely vytvorené pomocou XGBoost v R. SHAP hodnoty poskytujú podrobné vysvetlenie predikcií modelu na úrovni jednotlivých pozorovaní, čo pomáha v interpretácii “black box” modelov. Umožnujú pochopiť, ako dôležité sú premenné v rámci modelu. doParallel slúži na paralelné výpočty, a teda skrátenie času trénovania modelu. Posledná knižnica extrafont slúži na načítanie fontov inštalovaných na Windows-e.

Knižnice 
# import libraries
if (!require("pacman")) {
  install.packages("pacman")
}

pacman::p_load(
  tidyverse,
  sf, # data contain geometry
  tidymodels,
  vip,
  shapviz,
  doParallel,
  extrafont,
  vetiver,
  xgboost
)

loadfonts(device = "win")
options(scipen = 999)

unregister_dopar <- function() {
  env <- foreach:::.foreachGlobals
  rm(list = ls(name = env), pos = env)
}
Dáta 
# import data
original_conditions <- c(
  "Pôvodný stav", "Čiastočná rekonštrukcia", "Kompletná rekonštrukcia",
  "Novostavba", "Vo výstavbe", "Developerský projekt"
)
english_conditions <- c(
  "Original condition", "Partial reconstruction", "Complete reconstruction",
  "New building", "Under construction", "Development project"
)

original_type <- c("3 izbový byt", "1 izbový byt", "2 izbový byt", "4 izbový byt", "Garsónka", "5 a viac izbový byt", "Dvojgarsónka")
english_type <- c("3-room apartment", "1-room apartment", "2-room apartment", "4-room apartment", "Studio", "5 or more room apartment", "Two-room apartment")

apartments_analysis_data <- readRDS("data/apartments_final_data.rds") |>
  filter(!is.na(price)) |>
  mutate(
    coord = st_coordinates(st_centroid(geometry)),
    lon = coord[, 1],
    lat = coord[, 2],
    type = str_replace_na(recode(type, !!!setNames(original_type, english_type)),"Neznáme"),
    condition = str_replace_na(recode(condition, !!!setNames(original_conditions, english_conditions)),"Neznáme"),
    certificate = str_replace_na(str_replace(certificate, "none", "Nemá"), "Neznáme")
  ) |>
  select(-coord) # Optionally remove the original coordinates column

# remove geometry since we have coordinates now
apartments_analysis_data$geometry <- NULL

saveRDS(apartments_analysis_data, "data/apartments_data_App.RDS")

Tidymodels framework


Tidymodels je súbor knižníc v R, ktorý poskytuje jednotné a flexibilné rozhranie pre celý proces strojového učenia, od predspracovania dát cez ich analýzu až po modelovanie. Vytvorenie modelu v rámci frameworku tidymodels pozostáva z niekoľkých základných krokov:

  • Príprava a rozdelenie dát: Prvý krok je rozdelenie nášho “dátového budgetu”. Typicky sa dáta delia na trénovaciu (slúži na odhad parametrov modelu) a testovaciu (slúži na nezávislé zhodnotenie modelu) sadu pomocou funkcie initial_split(). Tento krok je základom pre overovanie modelu a zabránenie overfittingu. S využitím vfold_cv() alebo podobných funkcií vytvoríme schému krížovej validácie (rôzne verzie tréningových dát – tzv. “folds”), ktorá sa použije na evaluáciu modelu.
    Ako rozdelenie na tréningovú a testovaciu sadu, tak aj vytvorenie validačných schém umožnuje specifikovať premennú, ktorej rozdelenie ostane (približne) zachované (strata =).

Rozdelenie dátového budgetu. Zdroj: https://www.tidymodels.org/start/resampling/#data-split

  • Vytvorenie receptu (recipe): recipe() definuje predspracovanie dát, vrátane výberu premenných, transformácií, normalizácie, kódovania kategorických premenných a riešenia chýbajúcich hodnôt. Jednotlivé kroky sa pridávajú pomocou step_*() funkcií. Recipies zabezpečujú, že predspracovanie je konzistentné a reprodukovateľné.

  • Špecifikácia modelu: Model sa špecifikuje nezávisle od dát. Pomocou funkcií ako linear_reg(), rand_forest(), boost_tree() a iných definujeme typ modelu, mód (regresia, klasifikácia), engine(xgboost, lightgbm…) a jeho hlavné parametre bez toho, aby sme ich ihneď fitovali na dáta. Tento krok umožňuje flexibilitu v experimentovaní s rôznymi modelmi.

  • Nastavenie workflow: workflow() integruje recept a model do jednotného objektu. Workflow umožňuje efektívnejšie spracovanie, keďže spojíme predspracovanie dát a modelovanie do jednej operácie, čo zjednodušuje evaluáciu a porovnávanie modelov.

  • Výber a nastavenie hyperparametrov: Pomocou parameters() môžeme definovať a prispôsobiť rozsahy hyperparametrov pre tuning modelu. Tidymodels ponúka rôzne metódy pre vyhľadávanie optimálnych hodnôt, napr. tune_grid(), tune_bayes(), tune_race() a iné.

  • Cross-validácia a tuning modelu: Tuning hyperparametrov prebieha na trénovacej sade (resp. na jednotlivých “fold-och”) s cieľom nájsť najlepšiu kombináciu hyperparametrov. Najlepší model vyberieme pomocou funkcie select_best(), pričom špecifikujeme metriku, podľa ktorej model vyberáme.

  • Finalizácia a fitovanie modelu: Po vybraní najlepších hyperparametrov finalizujeme model pomocou finalize_model() a potom ho fitujete na trénovacie dáta s fit(). Tento krok produkuje finálny model pripravený na evaluáciu a predikcie.

  • Evaluácia modelu: Pomocou testovacej sady dát overíme výkonnosť modelu. Metriky ako RMSE, presnosť (accuracy), AUC a mnoho iných (výber by mal zodpovedať nášmu cieľu, čo platí najmä pri klasifikácii) poskytujú hodnotenie, ako dobre model predpovedá nevidené dáta.

Tieto kroky poskytujú ucelený prístup k vytváraniu, optimalizácii a evaluácii prediktívnych modelov. Tidymodels zabezpečuje konzistenciu a reprodukovateľnosť po celom procese.

Toto je opis základného procesu tvorby ML modelu pomocou tidymodels frameworku. Je možné samozrejme vytvoriť aj omnoho komplikovanejší proces s postupným tuningom hyperparametrov, tvorbou stacked modelu (meta-learner modelu) atď.

Tvorba modelu

Rozdelenie dát 
# split dataframes to train(80)/test(20)
set.seed(123)
apartments_train_split <- initial_split(apartments_analysis_data, prop = 0.7, strata = price)

apartments_train <- training(apartments_train_split)
apartments_test <- testing(apartments_train_split)

folds <- vfold_cv(apartments_train, v = 5, strata = price)
saveRDS(apartments_test, "data/test_set.RDS")
Spracovanie dát 
apartments_xgboost_recipe <- recipe(apartments_train, price ~ .) |>
  step_rm(name_nsi) |>
  step_string2factor(all_nominal_predictors()) |>
  step_impute_knn(all_nominal_predictors()) |>
  step_unknown(all_nominal_predictors()) |>
  step_dummy(all_nominal_predictors())
Špeficikácia modelu 
xgb_model <-
  boost_tree(
    trees = tune(), loss_reduction = tune(),
    tree_depth = tune(), min_n = tune(),
    mtry = tune(), sample_size = tune(),
    learn_rate = tune()
  ) |>
  set_mode("regression") |>
  set_engine("xgboost")
Workflow 
apartments_workflow <- workflow() |>
  add_recipe(apartments_xgboost_recipe) |>
  add_model(xgb_model)
Výber parametrov 
apartments_xgboost_params <- parameters(
  trees(), learn_rate(), loss_reduction(),
  tree_depth(), min_n(),
  sample_size = sample_prop(),
  finalize(mtry(), apartments_train)
)

apartments_xgboost_params <- apartments_xgboost_params |> update(trees = trees(c(300, 600)))


Pri tuningu použijeme paralelné spracovanie na x – 2 jadrách. Niektoré algoritmy umožnujú aj spracovanie pomocou GPU ak je dostupná kompatibilná grafická karta (bavíme sa o tvorbe modelu na lokálnom zariadení, samozrejme v produkcii je vhodnejšie využiť služby ako GCP, AWS a iné)

Tuning parametrov 
registerDoParallel(cores = detectCores() - 1)

xgboost_tune <-
  apartments_workflow |>
  tune_bayes(
    resamples = folds,
    param_info = apartments_xgboost_params,
    iter = 100,
    metrics = metric_set(rmse, mape),
    control = control_bayes(
      no_improve = 20,
      save_pred = T, verbose = F
    )
  )

unregister_dopar()
Finalizácia modelu 
apartments_best_model <- select_best(xgboost_tune, metric = "rmse")
apartments_final_model <- finalize_model(xgb_model, apartments_best_model)
apartments_workflow <- apartments_workflow |> update_model(apartments_final_model)
apartments_xgb_fit <- fit(apartments_workflow, data = apartments_train)

saveRDS(apartments_xgb_fit, "data/xgb_fit.RDS")
Výsledné metriky a testovacia predikcia 
apartments_final_res <- last_fit(apartments_workflow, split = apartments_train_split)
apartments_pred <-
  predict(apartments_xgb_fit, apartments_test) |>
  bind_cols(apartments_test)

Vyhodnotenie modelu

Z grafu je zrejmé, že výsledný model vcelku dobre predikuje ceny bytov s hodnotou do približne 300 tisíc EUR. Pri drahších bytoch je variabilita odchýlok od reálnych cien vyššia.

Porovnanie predikcie s reálnou cenou 1 
plot1 <-
  apartments_pred |>
  ggplot(aes(x = .pred, y = price)) +
  geom_point() +
  geom_smooth(method = "loess", color = "red") +
  scale_y_continuous(labels = function(x) format(x, big.mark = " ", scientific = FALSE), breaks = c(200000, 400000, 600000)) +
  scale_x_continuous(labels = function(x) format(x, big.mark = " ", scientific = FALSE), breaks = c(200000, 400000, 600000)) +
  labs(
    title = NULL,
    x = "Predikovaná cena",
    y = "Reálna Cena"
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "Courier New", size = 12)
  ) +
  coord_fixed()

plot1

Korelácia reálnych a predikovaných cien

Pri pohľade na rozloženie reálnych a predikovaných cien môžeme vidieť kde presne dochádza k najväčším odchýlkam. Model predikoval viac bytov v cene približne 100 až 170 tisíc EUR a približne 200 až 225 tisíc EUR. Naopak menej zastúpené sú najmä byty s cenou približne 250 až 300 tisíc EUR. Teraz je potrebné rozhodnutie, či sa treba vrátiť k procesu feature engineering, čiže manipulácii s premennými alebo pokračovať s model v stave v akom je.

Porovnanie predikcie s reálnou cenou 2 
plot2 <-
  apartments_pred |>
  select(predikcia = .pred, realita = price) |>
  gather(key, value) |>
  ggplot(aes(x = value, color = key)) +
  geom_density(alpha = .5) +
  labs(
    title = NULL,
    x = "Cena",
    y = "Hustota",
    color = ""
  ) +
  scale_y_continuous(labels = function(x) format(x, big.mark = " ", scientific = FALSE)) +
  scale_x_continuous(labels = function(x) format(x, big.mark = " ", scientific = FALSE)) +
  theme_minimal() +
  theme(
    text = element_text(family = "Courier New", size = 12),
    legend.position = "bottom"
  )

plot2

Rozloženie reálnych vs predikovaných cien

Model v priemere mierne nadhodnocuje (do 5 percent) byty s reálnou cenou medzi 100 až 170 tisíc EUR. Najvyššie nadhodnotenie (viac ako 5 percent) je pri cenách do 100 tisíc EUR. Naopak byty s cenou od 180 tisíc EUR do 260 tisíc eur mierne podhodnocuje. V prípade bytov od 260 tisíc EUR (približne 7 percent z inzerátov) sú už odchýlky nestabilné, trend je však smerom k výraznejšiemu podhodnocovaniu.

Odchýlky 
# Calculate the difference between prediction and actual value
plot3_data <- apartments_pred |>
  select(.pred, price) |> 
  mutate(diff = .pred - price,
         bin = floor(price/10000)*10000) |>
  group_by(bin) |>
  summarize(mean_diff = mean(diff)) |> 
  mutate(rel_diff = case_when(
    mean_diff/bin >= 0.05 ~ 0.05,
    mean_diff/bin >= 0 ~ 0,
    mean_diff/bin <= -0.05 ~ -0.05,
    mean_diff/bin < 0 ~ -0.02
    ))


# Plot with coloring based on mean difference
plot3 <- ggplot(plot3_data, aes(x = bin, y = mean_diff, color = factor(rel_diff))) +
  geom_point() +
  labs(
    title = NULL,
    x = "Cena",
    y = "Priemerná chyba"
  ) +
  scale_color_manual(values = c("#506BA0", "#55B9F5", "#90ee90","#037f51" )) +
  scale_y_continuous(labels = function(x) format(x, big.mark = " ", scientific = FALSE)) +
  scale_x_continuous(labels = function(x) format(x, big.mark = " ", scientific = FALSE)) +
  theme_minimal() +
  theme(
    text = element_text(family = "Courier New", size = 12),
    legend.position = "none")

plot3

Odchýlky predikcie od reality

Pri XGBoost je dôležitosť premenných založená na Gain, tento termín odkazuje na príspevok každej premennej k celkovému zlepšeniu modelu, ktoré je dosiahnuté vďaka rozdeleniam (splits) na konkrétnej premennej. Keď XGBoost vytvára stromy, každé rozdelenie v strome sa vyberá tak, aby maximalizovalo “zisk”, čo znamená, že sa snaží o čo najväčšie zlepšenie prediktívnej presnosti. Gain v tomto kontexte meria, o koľko sa zlepšuje predikcia, keď sa použije rozdelenie založené na danej premennej. Tento zisk je často vážený a sumarizovaný cez všetky stromy v modeli. Premenné, ktoré najviac prispievajú k zlepšeniu sú tie, ktoré model najviac využíva na dosiahnutie presnejších predikcií, a preto sú považované za dôležitejšie. Viac informácií môžete nájsť napríklad v tomto článku na medium.com.

Najdôležitejšou premennou je podľa modelu geografická dĺžka. Nasledujú rozloha bytu a lokalita Bratislava I. Až 4 z top 10 premenných súvisia s indexom bývania.

vip – Variable Importance Plots 
apartments_xgb_fit |>
  fit(data = apartments_train) |>
  pull_workflow_fit() |>
  vip(geom = "point", include_type = TRUE) + 
  theme_minimal() +
  theme(
    panel.grid.minor = element_blank(),
    text = element_text(family = "Courier New", size = 12)
  ) + 
  labs(
    y = "Dôležitosť",
    x = "Premenná"
  ) +
  theme(
    panel.grid.minor = element_blank(),
    text = element_text(family = "Courier New", size = 12)
  )

Dôležitosť premenných v modeli


SHAP (Shapley Additive exPlanations) hodnoty sú metóda používaná na vysvetlenie príspevku jednotlivých premenných k predikcii konkrétneho modelu strojového učenia. SHAP hodnoty poskytujú detailné vysvetlenie predikcie pre každý záznam (riadok dát) tým, že ukážu, ako každá premenná prispieva k výslednej predikcii, či už zvyšovaním alebo znížením predikovanej hodnoty. Tiež zaručujú konzistentnosť, čo znamená, že ak máme dve premenné a jedna konzistentne prispieva k predikcii viac než druhá, bude mať aj vyššiu SHAP hodnotu. Okrem poskytovania vysvetlení na úrovni jednotlivých predikcií (lokálne vysvetlenie) môžu byť SHAP hodnoty agregované na poskytnutie prehľadu o dôležitosti premenných v celom modeli (globálne vysvetlenie).

SHAP hodnoty nám poskytujú ďaľšie informácie k dôležitosti premenných. Napr. sa potvrdzuje západo-východný gradient, ktorý sme spomínali v závere tretej časti tejto série. Ak sa byt nachádza v okrese Bratislava I, môže mu to pridať na hodnote 45 až 95 tisíc EUR. Ak sa jedná o novostavbu cena môže stúpnuť o 15 až 50 tisíc EUR, kompletná rekonštrukcia pridáva na hodnote do 25 tisíc EUR, naopak pôvodný stav má negatívny efekt do približne -20 tisíc EUR.

SHAP 
fitted_data <- apartments_xgboost_recipe |>
  prep() |>
  bake(new_data = apartments_analysis_data) |>
  select(-price)

set.seed(123)
shp <- shapviz(extract_fit_engine(apartments_xgb_fit), X_pred = fitted_data |> as.matrix())

shap_plot <- sv_importance(shp, kind = "beeswarm")  +
  scale_x_continuous(labels = function(x) format(x, big.mark = " ", scientific = FALSE)) +
  scale_color_gradient(low = "blue", high = "red",
                      breaks = c(0,1),
                      labels = c("nízka", "vysoká"),
                      guide = guide_colorbar(barwidth = 12, barheight = 0.3)) + # Customize based on your exact needs
  labs(x = "SHAP hodnota (vplyv na predikovanú cenu)",
       y = "Premenná",
       color = "Hodnota premennej "
       ) +
  theme_minimal() +
  theme(
    panel.grid.minor = element_blank(),
    text = element_text(family = "Courier New", size = 10),
    legend.position = "bottom"
  )

Top 15 premenných podľa SHAP hodnôt

Okrem celkového príspevku jednotlivých premenných v rámci modelu ako celku (globálne vysvetlenie) sa vieme pozrieť aj na jednotlivé predikcie. Napr. keď sa pozrieme na inzerát na riadku 1804, môžeme vidieť, že lokalita v okrese Bratislava I mu pridala na hodnote 63 tisíc EUR, geografická dĺžka ďalších 33 tisíc EUR a rozloha 84 m2 20 tisíc EUR. Naopak to, že sa nejedná o novostavbu znižuje cenu o 12 tisíc EUR a pôvodný stav uberá ďalších 9 tisíc EUR.

Príklad SHAP hodnôt jednej predikcie

Uloženie modelu a dát pre aplikáciu

Aby sa model dal použiť v Shiny aplikácii, je nutné ho uložiť vo vhodnom formáte. Natrénovaný model transformujeme na vetiver model. Tento objekt obsahuje všetko, čo treba, aby sa dal použiť v novom prostredí. Ďalej ukladáme engine pomocou knižnice xgboost. Budeme ho potrebovať pri SHAP grafe.