Machine learning v jazyku R – Odhad cien bytov

Ú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.

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.

# 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)
}

# 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 =).

• 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

# 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")

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

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")

apartments_workflow <- workflow() |>
  add_recipe(apartments_xgboost_recipe) |>
  add_model(xgb_model)

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é)

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

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")

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.

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

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.

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

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.

# 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

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.

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

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.

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"
  )

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 m² 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.

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.