Themen und Lernziele

Nachdem wir in der Einführung versuchen das Thema im Allgemeinen zu verstehen und uns mit statistischen Konzepten überlegen, warum Ensembles oft sehr gut funktionieren, werden wir in den drei nächsten Kapiteln drei populäre Ensemble Methoden anschauen: Bagging, Random Forests und Boosting.

Warum Ensembles?

Doch warum funktionieren Ensembles so gut? Hier hilft eine Analogie zur klassischen Statistik (ISLR, p. 340). Wir nehmen an, wir haben eine Stichprobe von \(n\) Beobachtungen (z.B. Hauspreise).

• Die Standardabweichung \(\sigma\) besagt, wie stark eine einzelne Beobachtung im Schnitt vom Stichprobenmittelwert abweicht.
• Der Standardfehler des Mittelwerts \(\frac{\sigma}{\sqrt{n}}\) besagt, wie stark der Stichprobenmittelwert im Schnitt vom wahren (unbekannten) Populationsmittelwert abweicht.

Wir sehen, dass sich die Varianz durch das Mitteln um den Faktor \(1/\sqrt{n}\) reduziert.

Bei Ensemble Methoden passiert etwas ähnliches: indem wir die Vorhersagen von vielen individuellen Modellen mitteln, reduziert sich die Varianz. Basierend auf dem Bias-Variance Tradeoff können wir ableiten, dass ein Ensemble Modell mit geringerer Varianz auch zu besseren Vorhersagen führen sollte.

Wichtig: die Formel für den Standardfehler des Mittelwerts in der klassischen Statistik gilt nur, wenn die \(n\) Beobachtungen unabhängig sind. Auch dieses Resultat lässt sich auf Ensemble Methoden übertragen: je unabhängiger die individuellen Modelle voneinander sind, desto besser funktioniert ein Ensemble.

Voting Models

In vielen Projekten trainieren wir mehrere Modelle und entscheiden uns dann für das beste Modell. Doch warum nicht alle trainierten Modelle verwenden? Genau das ist die Idee von sogenannten Voting Models. Die Idee könnte simpler nicht sein: wir lassen jedes unserer Modelle eine Vorhersage machen (harte Vorhersage oder Wahrscheinlichkeit) und kombinieren alle individuellen Vorhersagen zu einer aggregierten Vorhersage. Das Prinzip des Voting Models ist in folgender Abbildung anhand des Spam Filters illstriert:

Classification Trees auf dem Heart Datensatz

Hier wollen wir nochmals einen Classification Tree auf dem Heart Datensatz rechnen und die Resultate dann später mit den Ensemble Methoden vergleichen.

Wir laden als erstes diverse R-Packages und danach die Daten. Im Vergleich zum letzten Mal haben wir die Beobachtungen mit fehlenden Werten in mindestens einer Spalte entfernt, weil die Modelle, die wir später anschauen werden, nicht gut mit fehlenden Werten umgehen können. Darum enthält der heart Dataframe nur noch 297 Zeilen.

# Laden der Libraries
library(rpart.plot)
library(tidyverse)
library(tidymodels)
library(ranger)
library(randomForest)

# Tidy Workspace
rm(list = ls())

# Working Directory
setwd("...")

# Daten importieren
heart <- read_rds("heart.rds")

Wir machen hier nun einen Train-Test Split, weil wir dann später alle Modelle auf einem Testdatensatz evaluieren wollen. Wir bilden auch schon gleich die 5 Folds für die Cross-Validation.

# ------------------------------------------------------
# 1. Vorarbeiten

# Ist der Datensatz balanced?
heart %>% 
  count(hd) %>% 
  mutate(prop = n / sum(n))

# Seed, um den Train-Test Split reproduzierbar zu machen
set.seed(123)

# Split (75% Training, 25% Test)
train_test_split <- initial_split(heart, prop = 3/4, strata = hd)

train <- training(train_test_split)
test  <- testing(train_test_split)

# Seed, um 5-Fold Cross-Validation reproduzierbar zu machen
set.seed(456)

# 5-Fold CV
folds <- vfold_cv(train, v = 5, repeats = 1, strata = hd)

Als nächstes nun der Code, um einen Classification Tree zu fitten. Der Code ist grösstenteils identisch zum Code vom letzten Mal:

# ------------------------------------------------------
# 2. Classification Tree

# Spezifikation des Classification Trees
dt_mod <- 
  decision_tree(tree_depth = tune(), min_n = tune(), cost_complexity = 0) %>%
  set_mode("classification") %>% 
  set_engine("rpart")

# Workflow
dt_workflow <- 
  workflow() %>% 
  add_model(dt_mod) %>% 
  add_formula(hd ~ .)

# Tuning Grid
dt_grid <- expand.grid(tree_depth = seq(1, 20, by = 2), min_n = seq(10, 30, 5))

# Tuning / Model Fitting
dt_res <- 
  dt_workflow %>% 
  tune_grid(
    resamples = folds,
    grid = dt_grid,
    control = control_grid(save_pred = TRUE),
    metrics = metric_set(roc_auc, accuracy))

# Wir sortieren die Hyperparameter Spezifikationen nach ROC AUC
dt_res %>% 
  show_best(metric = "roc_auc", n = 10) %>% 
  arrange(desc(mean))

# Wir wählen das beste Modell
dt_best <- 
  dt_res %>% 
  select_best(metric = "roc_auc")

# Resultate vorbereiten für ROC Kurve
dt_auc <- 
  dt_res %>% 
  collect_predictions(parameters = dt_best) %>% 
  roc_curve(hd, .pred_yes, event_level = "second") %>% 
  mutate(model = "Classification Tree")

# ROC Kurve
autoplot(dt_auc)

# Accuracy
dt_res %>% 
  collect_predictions(parameters = dt_best) %>% 
  accuracy(hd, .pred_class)

# Last Fit Spezifikation (optimale Hyperparameter Werte)
last_dt_mod <-
  decision_tree(tree_depth = 5, min_n = 20, cost_complexity = 0) %>% 
  set_mode("classification") %>% 
  set_engine("rpart")

# Workflow anpassen (optimales Modell)
last_dt_workflow <- 
  dt_workflow %>% 
  update_model(last_dt_mod)

# Last Fit (auf ganzem Trainingsdatensatz)
last_dt_fit <-
  last_dt_workflow %>% 
  fit(data = train)

# Last Fit
last_dt_fit

# Plot des finalen Baums
last_dt_fit %>%
  extract_fit_engine() %>%
  rpart.plot()

Wenn wir obigen Code in R ausführen, dann kriegen wir mithilfe der Cross-Validation bereits eine Schätzung der Accuracy. Die Accuracy basierend auf den CV-Resultaten beträgt 79.3%. Dementsprechend beträgt die Error Rate 20.7%. Es ist hier akzeptabel die Accuracy zu verwenden, da der Datensatz mehr oder weniger ausgeglichen ist.

Weiterführende Ressourcen

• Kapitel 8 in ISLR
• Kapitel 7 in HOML
• Dokumentation baguette
• Random Forest Dokumentation tidymodels
• Bagging with tidymodels and #TidyTuesday astronaut missions
• Tuning random forest hyperparameters with #TidyTuesday trees data
• A Gentle Introduction to tidymodels
• Gradient Boosting Machines
• XGBoost Demos
• Tune XGBoost with tidymodels and #TidyTuesday beach volleyball

Bootstrap

Die erste wichtige Frage ist, wie wir aus dem Trainingsdatensatz verschiedene Subsets generieren. Dazu wird beim Bagging eine altbekannte Technik aus der Statistik angewendet, nämlich der Bootstrap. Die beiden folgenden Abbildungen zeigen, wie der Bootstrap in der klassischen Statistik angewendet werden.

Der zentrale Punkt, den es hier zu verstehen gilt, ist wie das Boostrap Sampling funktioniert. Wir ziehen \(n\) Beobachtungen aus einer Stichprobe der Grösse \(n\). Damit jedes Subset (d.h. jedes Bootstrap Sample) unterschiedlich ist, müssen wir ziehen mit Zurücklegen. Dadurch werden gewisse Beobachtungen mehr als einmal in unserem Bootstrap Sample landen und gewisse Beobachtungen gar nicht. Nur so kriegen wir unterschiedliche Bootstrap Samples.

In obigem Beispiel können wir den Standardfehler des Mittelwerts \(s_{\bar{x}}\) selbstverständlich analytisch berechnen, nämlich mithilfe der Formel \(s_{\bar{x}}=\frac{\hat{\sigma}}{\sqrt{n}}=\frac{316.47}{\sqrt{10}}\approx 100\). Manchmal gibt es aber keine analytische Formel und dann ist die Boostrap Technik enorm hilfreich. Die nächste Abbildung zeigt, wie wir den Standardfehler des Mittelwerts mithilfe des Boostraps rechnen könnten:

Die Standardabweichung über die 1000 Mittelwerte, die wir basierend auf den Bootstrap Samples rechnen, gibt uns eine Schätzung für den Standardfehler des Mittelwerts. Unsere Bootstrap Schätzung (95) ist nicht weit entfernt von der analytischen Lösung (100)!

Wie funktioniert Bagging?

Wir generieren also mit der Bootstrap Technik eine grosse Anzahl \(B\) Subsets aus dem Trainingsdatensatz und rechnen für jedes Subset einen Decision Tree. Beim Bagging beschränken wir die Bäume nicht (keine max. Tiefe oder ähnliches). Jeder einzelne Baum wird darum viel Varianz, aber relativ wenig Bias haben.

Die Idee des Ensembles ist nun, dass wir die \(B\) Vorhersagen, die wir von all den individuellen Decision Trees kriegen, mitteln. Im Fall des Regressionsproblems rechnen wir ganz einfach den Durchschnitt über die \(B\) Vorhersagen,

\[ \hat{f}(\mathbf{x}_i) = \frac{1}{B} \sum_{b=1}^B \hat{f}^b(\mathbf{x}_i) \] wobei \(\hat{f}^b(\mathbf{x}_i)\) die Vorhersage des Baums im \(b\)-ten Subset bezeichnet.

Im Fall des Klassifikationsproblems können wir zwischen einem Soft Vote oder einem Hard Vote wählen:

• Beim Soft Vote rechnen wir den Durchschnitt über die vorhergesagten Wahrscheinlichkeiten der \(B\) Bäume.
• Beim Hard Vote entscheiden wir uns für die “harte” Vorhersage, die unter den \(B\) Bäumen am häufigsten vorkommt.

Folgende Abbildung zeigt das Vorgehen schematisch anhand des housingrents Datensatzes, den Sie im Decision Tree Tutorial bereits kennen gelernt haben. Es handelt sich hierbei um ein Regressionsproblem:

Wir aggregieren hier \(B=3\) Bäume zu einem Bagging Ensemble.

Out-of-Bag

Das Bootstrap Sampling führt dazu, dass in jedem Subset im Schnitt nur rund 2/3 der Beobachtungen im Trainingsdatensatz sind, d.h. mindestens einmal gezogen werden. Die nicht gezogenen Beobachtungen nennen wir die Out-of-Bag (OOB) Beobachtungen.

Warum ist OOB wichtig? Weil jede Beobachtung im Trainingsdatensatz in rund 1/3 aller Subsets OOB sein wird. Das können wir für eine Art Cross-Validation nutzen. Wir schauen für jede Beobachtung, in welchen Subsets sie OOB war. Dann verwenden wir die Bäume aller dieser Subsets, um für die OOB-Beobachtung Vorhersagen zu rechnen. Die Vorhersagen werden gemittelt. Wenn wir dies für jede Beobachtung im Trainingsdatensatz machen, dann kriegen wir für jede Beobachtung eine Vorhersage von Modellen, für welche diese Beobachtung nicht im Training verwendet wurde.

Für grosse Trainingsdatensätze kann das ein sehr wertvolles Vorgehen sein, weil K-Fold Cross-Validation in diesen Fällen sehr “teuer” ist (R läuft sehr lange!).

Variable Importance

Wir bewegen uns nun immer stärker in Richtung Black-Box Modelle. Darum wird das Thema Variable Importance (VIP) nun wichtiger. Das heisst, wir müssen uns nun verstärkt überlegen, wie wir aus dem Modell lernen können, welche Input-Variablen für die Vorhersagen wichtig sind.

Beim linearen und logistischen Regressionsmodell war das relativ einfach, da man die Wichtigkeit einer Input-Variable sehr gut anhand des geschätzten Koeffizienten für die Variable abschätzen kann. Auch bei Decision Trees ist es simpel: Variablen, welche in Splits weiter oben im Baum verwendet werden, sind tendenziell wichtiger als Variablen, die in Splits weiter unten vorkommen und garantiert wichtiger als Variablen, die im Baum nicht verwendet werden.

Die gute Nachricht ist, dass es auch für Ensembles gute und konzeptionell einfache Methoden gibt. Bagging ist ja bekanntlich ein Modell, das über eine grosse Anzahl \(B\) Decision Trees aggregiert. Wenn wir beispielsweise die Wichtigkeit für die Input-Variable \(x_{i1}\) messen wollen, dann schauen wir in jedem der \(B\) Bäume, wie gross die Reduktion in SQR (Regression) bzw. Gini Index (Klassifikation) ist durch einen Split, in dem die Variable \(x_{i1}\) verwendet wurde. Die VIP der Variable \(x_{i1}\) ist dann einfach der Durchschnitt über die \(B\) Reduktionen (wenn die Variable \(x_{i1}\) in einem der Bäume nicht verwendet wird, dann beträgt die Reduktion für diesen Baum 0).

Häufig rechnet man VIP gleich für alle Input-Variablen und normalisiert sie, so dass VIP ein Maximum von 1 annehmen kann. Das erreichen wir, indem wir alle (unnormalisierten) VIP Werte durch den maximalen VIP Wert teilen.

Bagging auf dem Heart Datensatz

Wie oben schon etwas angetönt, ist es beim Bagging nicht zwingend nötig irgendwelche Hyperparameter zu tunen. Die Anzahl Bäume \(B\) sollte einfach hoch genug gewählt werden und die individuellen Bäume müssen nicht beschränkt werden. Aus diesem Grund werden wir hier nicht tidymodels verwenden (wir sparen uns hier den ganzen Overhead von tidymodels). Stattdessen verwenden wir das R Package randomForest. Sie werden später sehen, dass wir dasselbe Package auch für die Random Forest Methode verwenden werden (wie es der Name natürlich bereits verrät).

Wir übergeben der randomForest() Funktion als erstes die Modellformel hd ~ . und den Trainingsdatensatz. Das Argument mtry wird vor allem bei den Random Forest Modellen wichtig sein. Es gibt der Funktion vor, wie viele der Input-Variablen an jedem Split eines individuellen Baums als mögliche Splitvariablen zur Verfügung stehen. Hier, beim Bagging, sollen bei jedem Split alle Input-Variablen zur Verfügung stehen. Darum setzen wir das Argument auf ncol(train) - 1 (13), die Anzahl Input-Variablen. Das Argument ntree = 10000 gibt der Funktion vor, dass wir \(B=10'000\) Bäume schätzen wollen. Zum Schluss sagen wir noch, dass wir die Variable Importance messen wollen.

# ------------------------------------------------------
# 3. Bagging

# Seed für Reproduzierbarkeit (Bagging basiert stark auf SAMPLING)
set.seed(44)

# Fitting eines Bagging Modells
bag_fit <- 
  randomForest(hd ~ ., 
               data = train, 
               mtry = ncol(train) - 1, 
               ntree = 10000, 
               importance = TRUE)

# Output des Modells
bag_fit

Call:
 randomForest(formula = hd ~ ., data = train, mtry = ncol(train) - 1, ntree = 10000, importance = TRUE) 
               Type of random forest: classification
                     Number of trees: 10000
No. of variables tried at each split: 13

        OOB estimate of  error rate: 22.97%
Confusion matrix:
    no yes class.error
no  97  23   0.1916667
yes 28  74   0.2745098

Sie sehen oben eine Zusammenfassung des Fit Objekts. Insbesondere interessant ist die Out-of-Bag (OOB) Error Rate, die hier 22.97% beträgt und damit leicht höher ist als die CV Error Rate unseres Classification Trees (20.7%).

Selbstverständlich könnten wir das gefittete Modell nun bereits verwenden, um harte oder weiche Vorhersagen zu machen:

# Vorhersagen machen
predict(bag_fit, newdata = train[1, ], type = "prob")
predict(bag_fit, newdata = train[1, ], type = "class")

## Soft (type = "prob")
      no    yes
1 0.9845 0.0155
attr(,"class")
[1] "matrix" "array"  "votes" 

## Hard (type = "class")
 1 
no 
Levels: no yes

Wichtig: die predict() Funktion in obigem Code Block ist aus Base R. Bei der predict() Funktion aus tidymodels würde das Argument newdata anders heissen, nämlich new_data.

Zum Schluss können wir uns noch die Variable Importance anschauen:

# Variable Importance
randomForest::importance(bag_fit)

                no        yes MeanDecreaseAccuracy MeanDecreaseGini
age      17.100347   3.715216            16.177084        7.8986270
sex      31.081711   6.019096            29.483147        1.5685783
cp       40.770308  48.529442            61.272832       13.5563437
trestbps 14.167571   3.402533            13.105737        8.6907289
chol      4.071502 -15.546766            -6.707937        7.6470803
fbs      -2.062854  -6.660895            -5.817300        0.3817802
restecg  -5.875380  15.080576             7.003098        1.7126149
thalach  38.111493  19.476168            42.650182       12.2312568
exang    17.678602  10.001055            20.078778        1.8729007
oldpeak  82.348719  73.489513           107.656615       15.3368627
slope    13.423276  24.646575            27.516395        3.0001435
ca       94.181277  43.425076            95.753121       12.4442411
thal     83.619876  73.070677           102.413468       23.4377673

Wir sehen, dass die Variable thal klar zur stärksten durchschnittlichen Reduktion im Gini Index führt! Wir können diese Werte natürlich auch grafisch darstellen:

# Automatischer Plot Variable Importance
varImpPlot(bag_fit, main = "Variable Importance Heart")

Wie funktionieren Random Forests?

Um mehr Diversität (i.e., weniger Korrelation) unter den Bäumen herzustellen, macht die Random Forest Methode einen kleinen, aber cleveren Trick: für jeden individuellen Baum und jeden Split stehen nur \(m\) zufällig ausgewählte Input-Variablen zur Verfügung (wobei \(m\) kleiner ist als die Anzahl Input-Variablen \(p\)). So ergeben sich vielfältigere Bäume und die Reduktion der Varianz ist grösser. Etwas konkreter: die dominante Variable econage steht so nicht für jeden ersten Split zur Verfügung und so zwingen wir die Bäume zu mehr Diversität.

Typischerweise tunen wir \(m\) via Cross-Validation. In der Praxis hat sich \(m=\sqrt{p}\) als guter Wert herausgestellt. Wenn wir also kein Hyperparameter-Tuning machen wollen, dann kann \(m\) manuell auf \(\sqrt{p}\) gesetzt werden. Die Anzahl Bäume \(B\) kann problemlos auf einen grossen Wert gesetzt werden (z.B. \(B=10'000\)). Hierbei besteht keine Gefahr für Overfitting. Die Tatsache, dass ein Random Forest in vielen Fällen auch ohne Hyperparameter Tuning gute Resultate liefert, ist einer der Hauptgründe, warum dieses Modell in der Praxis so beliebt ist.

Ein Nachteil von Random Forests ist, dass die Trainingsphase länger dauern kann, insbesondere wenn wir uns entscheiden, gewisse Hyperparameter doch zu tunen. Wenn wir z.B. \(B=1000\) setzen, dann heisst das, dass wir für jeden Hyperparameter Wert 1000 Decision Trees schätzen müssen. Doch auch hier gibt es eine Lösung, denn die einzelnen Bäume werden unabhängig voneinander trainiert und darum ist es sehr einfach den Rechenprozess zu parallelisieren (auf mehrere Cores Ihres Computers).

Random Forests auf dem Heart Datensatz

Wir schauen uns hier zuerst an, wie man Random Forests ohne tidymodels rechnen kann. Dazu verwenden wir wie bereits beim Bagging die randomForest() Funktion. Der einzige Unterschied im Vergleich zum Bagging ist, dass wir nun für das Argument mtry einen Wert kleiner als die Anzahl Input-Variablen setzten. Wie oben erwähnt gilt die Faustregel, dass \(m=\sqrt{p}=\sqrt{13}=3.61\). Darum setzten wir mtry = 4. Nun stehen bei jedem Split in jedem individuellen Baum nur vier zufällig ausgewählte Input-Variablen zur Verfügung. Wie gesagt: ein kleiner, aber feiner Trick!

# ------------------------------------------------------
# 4.1 Random Forests (klassisch)

# Seed für Reproduzierbarkeit
set.seed(22)

# Fitting eines Random Forest Modells
rf_fit <- 
  randomForest(hd ~ ., 
               data = train, 
               mtry = 4, 
               ntree = 10000, 
               importance = TRUE)

# Output des Modells
rf_fit

Call:
 randomForest(formula = hd ~ ., data = train, mtry = 4, ntree = 10000, importance = TRUE) 
               Type of random forest: classification
                     Number of trees: 10000
No. of variables tried at each split: 4

        OOB estimate of  error rate: 19.37%
Confusion matrix:
     no yes class.error
no  103  17   0.1416667
yes  26  76   0.2549020

Wir sehen, dass die OOB Error Rate nun nur 19.37% beträgt. Wir dürfen also erwarten, dass der Random Forest bessere Vorhersagen machen wird als Bagging und der Classification Tree.

Nun fitten wir den Random Forest mit tidymodels und der ranger Engine. Schauen Sie sich folgenden Code genau an. Viele Elemente kennen Sie bereits. Was ist neu? Wir verwenden das R-Package doParallel, um das Tuning zu parallelisieren. Das heisst, jede Hyperparameter Kombination wird auf einem anderen Core Ihres Computers gerechnet. So verläuft der Trainingsprozess relativ schnell! Mit grid = 30 lassen wir einfach 30 Hyperparamter Kombinationen testen; wir geben sie nicht explizit vor.

# ------------------------------------------------------
# 4.2 Random Forests (tidymodels + ranger)

# Spezifikation des RFs
rf_mod <- 
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) %>% 
  set_engine("ranger") %>% 
  set_mode("classification")

# Workflow
rf_workflow <- 
  workflow() %>% 
  add_model(rf_mod) %>% 
  add_formula(hd ~ .)

# Parallel Processing
doParallel::registerDoParallel()

# Seed für Reproduzierbarkeit
set.seed(22)

# Model Fitting / Tuning
rf_res <- 
  rf_workflow %>% 
  tune_grid(folds,
            grid = 30,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))

# Wir sortieren die Hyperparameter Spezifikationen nach ROC AUC
rf_res %>% 
  show_best(metric = "roc_auc", n = 10) %>% 
  arrange(desc(mean))

# A tibble: 10 × 8
    mtry min_n .metric .estimator  mean     n std_err .config              
   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
 1     1     7 roc_auc binary     0.909     5  0.0274 Preprocessor1_Model22
 2     2     4 roc_auc binary     0.903     5  0.0278 Preprocessor1_Model25
 3     2    31 roc_auc binary     0.898     5  0.0302 Preprocessor1_Model02
 4     4    10 roc_auc binary     0.895     5  0.0264 Preprocessor1_Model16
 5     3    34 roc_auc binary     0.894     5  0.0296 Preprocessor1_Model09
 6     3    25 roc_auc binary     0.894     5  0.0301 Preprocessor1_Model23
 7     3    37 roc_auc binary     0.893     5  0.0293 Preprocessor1_Model24
 8     4    24 roc_auc binary     0.893     5  0.0280 Preprocessor1_Model27
 9     7    13 roc_auc binary     0.891     5  0.0245 Preprocessor1_Model04
10     5    18 roc_auc binary     0.889     5  0.0288 Preprocessor1_Model01

Interessanterweise ist für mtry der Wert 1 optimal. Das heisst für jeden Split steht nur eine zufällig ausgewählte Input-Variable zur Verfügung. Wir können die Tuning Resultate auch grafisch darstellen:

# Plot Tuning Resultate
autoplot(rf_res)

Nun speichern wir das beste Modell ab und vergleichen die ROC Kurven des Classification Trees und des Random Forests:

# Bestes Modell
rf_best <- 
  rf_res %>% 
  select_best(metric = "roc_auc")

# Resultate vorbereiten für ROC Kurve
rf_auc <- 
  rf_res %>% 
  collect_predictions(parameters = rf_best) %>% 
  roc_curve(hd, .pred_yes, event_level = "second") %>% 
  mutate(model = "Random Forest")

# ROC Kurven vergleichen
bind_rows(rf_auc, dt_auc) %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(linewidth = 1.5, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  scale_color_viridis_d(option = "plasma", end = .6) +
  theme_bw()

Der Random Forest ist doch einigermassen klar besser. Cool! Nun erstellen wir wie gewohnt den Last Fit auf dem ganzen Trainingsdatensatz und mit optimalen Hyperparameterwerten. Mit importance = "impurity" sagen wir der Funktion, dass wir die Variable Importance Werte auch gleich mitrechnen wollen.

# Last Fit Spezifikation (optimale Hyperparameter Werte)
last_rf_mod <-
  rand_forest(mtry = 1, min_n = 7, trees = 1000) %>% 
  set_mode("classification") %>% 
  set_engine("ranger", importance = "impurity")

# Workflow anpassen (optimales Modell)
last_rf_workflow <- 
  rf_workflow %>% 
  update_model(last_rf_mod)

# Last Fit (auf ganzem Trainingsdatensatz)
last_rf_fit <-
  last_rf_workflow %>% 
  fit(data = train)

Schauen wir uns kurz das gefittete Objekt an:

# Modell Output
last_rf_fit %>% 
  extract_fit_parsnip()

parsnip model object

Ranger result

Call:
 ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~1, x), 
 num.trees = ~1000, min.node.size = min_rows(~7, x), importance = ~"impurity", 
 num.threads = 1, verbose = FALSE, seed = sample.int(10^5, 1), probability = TRUE) 

Type:                             Probability estimation 
Number of trees:                  1000 
Sample size:                      222 
Number of independent variables:  13 
Mtry:                             1 
Target node size:                 7 
Variable importance mode:         impurity 
Splitrule:                        gini 
OOB prediction error (Brier s.):  0.1393055

Auch hier sehen wir den OOB Fehler. Er liegt bei knapp 14%. Hier handelt es sich jedoch um das sogenannte Brier Score, weshalb dieser Wert nicht direkt mit den vorherigen OOB Werten verglichen werden kann. Last but not least, können wir die Variable Importance Werte für die Top-20 Input-Variablen plotten mit der Funktion vip().

# Variable Importance Plot
last_rf_fit %>% 
  extract_fit_parsnip() %>% 
  vip(num_features = 20)

Die Variable thal ist in allen Modellen die klar wichtigste Variable.

Wie funktioniert Gradient Boosting?

Wir schauen uns hier der Einfachheit halber an, wie der Gradient Boosting Algorithmus für das Regressionsproblem funktioniert. Der Algorithmus für das Klassifikationsproblem ist ähnlich, aber technisch etwas komplexer.

Aus obigem Algorithmus ist ersichtlich, dass Gradient Boosting drei wichtige Hyperparameter hat:

• Die Anzahl Bäume \(B\), die in das Ensemble einfliessen. Im Gegensatz zu Bagging und Random Forest kann es beim Boosting zu Overfitting kommen, wenn wir eine zu grosse Anzahl Bäume \(B\) wählen, denn jeder Baum versucht ja die Fehler aller vorangehenden Bäume zu verbesseren oder korrigieren. Irgendwann kommen wir so in den Bereich, wo der Algorithmus nur noch den unsystematischen Teil der Information (Noise) zu fitten versucht.
• Der Hyperparameter \(d\) bezeichnet die Anzahl Splits, die wir in den einzelnen Bäumen zulassen. Erstaunlich oft funktionieren Bäume mit einem Split, also \(d=1\), gut in der Praxis. In diesem Fall handelt es sich um sogenannte Decision Stumps (dt. Entscheidungsstümmel). In diesem Fall fitten wir effektiv ein additives Modell (warum?).
• Der dritte wichtige Hyperparameter ist der Shrinkage Hyperparameter \(\lambda \leq 1\). Wir wollen damit sicherstellen, dass das Boosting Ensemble nicht zu schnell lernt. Es hat sich in der Praxis herausgestellt, dass langsam lernende Boosting Ensembles zu einem besseren Modell führen. Der Shrinkage Hyperparameter kann auch als eine Form von Regularisierung interpretiert werden.

Wir werden später im Deep Learning den Gradient Descent Algorithmus kennen lernen. In der folgenden optionalen Box lernen Sie, inwiefern Gradient Boosting mit Gradienten und Ableitungen zu tun hat.

Gradient Boosting auf dem Heart Datensatz

Wir fitten hier den Gradient Boosting Algorithmus in einem ersten Schritt mit dem traditionellen gbm Package. Dieses Package erwartet, dass die kategorische Zielvariable 0-1 kodiert ist, weshalb ich in einem ersten Schritt eine Umkodierung vornehme.

Danach fitten wir mit der Funktion gbm() das Modell. In der Modellformel entferne ich die ursprüngliche Zielvariable hd. Mit distribution = "bernoulli" sagen wir der Funktion, dass es sich um ein binäres Klassifikationsproblem handelt. Die drei wichtigen Hyperparameter heissen hier etwas anders:

• \(B\) → n.trees
• \(d\) → interaction.depth
• \(\lambda\) → shrinkage

Sie sehen, dass wir via dem Argument cv.folds = 5 auch eine interne Cross-Validation betreiben. Wir wollen damit herausfinden, wie viele Trees nötig sind und ob es allenfalls weniger Bäume als n.trees = 10000 braucht.

# ------------------------------------------------------
# 5.1 Gradient Boosting (gbm)

# Seed für Reproduzierbarkeit
set.seed(15)

# Wir müssen Zielvariable in 0/1 umkodieren für gbm Package
train$y <- as.integer(train$hd) - 1

# Fitting eines GBM Modells
gbm_fit <- gbm(
  formula = y ~ . - hd,
  distribution = "bernoulli",
  data = train,
  n.trees = 10000,
  interaction.depth = 1,
  shrinkage = 0.001,
  cv.folds = 5,
  n.cores = NULL,
  verbose = FALSE
)

# Output
print(gbm_fit)

gbm(formula = y ~ . - hd, distribution = "bernoulli", data = train, 
    n.trees = 10000, interaction.depth = 1, shrinkage = 0.001, 
    cv.folds = 5, verbose = FALSE, n.cores = NULL)
A gradient boosted model with bernoulli loss function.
10000 iterations were performed.
The best cross-validation iteration was 7362.
There were 13 predictors of which 13 had non-zero influence.

Ha! Tatsächlich ist die optimale Anzahl Bäume (also das optimale \(B\)) tiefer als 10’000 und liegt bei 7362 Bäumen. Der Cross-Validation Fehler ist bei dieser optimalen Anzahl Bäumen am tiefsten. Dieses Resultat kann man sich auch noch grafisch ausgeben lassen:

# Kostenfunktion für steigende Anzahl Trees
gbm.perf(gbm_fit, method = "cv")

Sie sehen, dass wir nach der blau gestrichelten Vertikalen in den Overfitting Bereich kommen: die Kosten auf dem Trainingsdatensatz (in schwarz) gehen weiter runter, aber die Kosten auf den CV Folds beginnen wieder leicht zu steigen (in anderen Worten: die Varianz des Modells beginnt zu steigen).

Wir fitten nun also das optimale Gradient Boosting Modell mit der optimalen Anzahl Bäumen. Die anderen Hyperparameter tunen wir hier der Einfachheit halber nicht. Mithilfe des Packages vip lassen wir uns die Wichtigkeit der verschiedenen Variablen visualisieren.

# Finales GBM Modell
gbm_final_fit <- gbm(
  formula = y ~ . - hd,
  distribution = "bernoulli",
  data = train,
  n.trees = 7362,
  interaction.depth = 1,
  shrinkage = 0.001,
  n.cores = NULL,
  verbose = FALSE
)

# Variable Importance
vip(gbm_final_fit)

Wie bereits vermutet ist auch hier die Variable thal am wichtigsten!

Nun werden wir uns anschauen, wie wir das Gradient Boosting Modell mit tidymodels und der xgboost Library fitten und tunen können. xgboost steht für Extreme Gradient Boosting und ist eine optimierte und enorm flexible Implementation des Gradient Boosting Algorithmus.

xgboost hat nicht nur die drei oben erwähnten Hyperparamter, sondern viele weitere, welche das Modell schneller und besser machen können. Eine Übersicht über die Hyperparameter finden Sie hier.

Im unten stehenden Code entfernen wir als erstes die 0/1 kodierte Zielvariable, die wir für die gbm Implementation benötigten. Danach spezifizieren wir in typischer tidymodels Manier das Modell und den Workflow. Der Einfachheit halber lassen wir nur die drei bekannten Hyperparameter trees (\(B\)), tree_depth (\(d\)) und learn_rate (\(\lambda\)) tunen. Die weiteren Hyperparameter belassen wir auf den Default Werten. Mit nthread = 8 erlauben wir dem Modell ein sogenanntes Parallel Processing (eines der tollen Features der xgboost Implementation).

In der tune_grid() Funktion setzen wir grid = 50 und spezifizieren keinen eigenen Grid mit möglichen Hyperparameterkombinationen. Dadurch erstellt tidymodels automatisch ein Space-filling Design. Doch was bedeutet das? Aus allen möglichen Kombinationen von Werten für trees, tree_depth und learn_rate wählt tidymodels automatisch 50 Kombinationen, die Sinn machen und die den Raum (Space), welcher durch die drei Hyperparameter aufgespannt wird, möglichst gut ausfüllen (Space-filling).

# ------------------------------------------------------
# 5.2 Gradient Boosting (xgboost)

# Wir entfernen y wieder aus Trainingsdatensatz
train <- train[ ,-ncol(train)]

# Spezifikation des XGBoost Modells
xgb_mod <- 
  boost_tree(trees = tune(), tree_depth = tune(), learn_rate = tune()) %>% 
  set_engine("xgboost", event_level = "second", nthread = 8) %>% 
  set_mode("classification")

# Workflow
xgb_workflow <- 
  workflow() %>% 
  add_model(xgb_mod) %>% 
  add_formula(hd ~ .)

# Seed für Reproduzierbarkeit
set.seed(42)

# Model Fitting / Tuning
# So genanntes "Space-filling design": es werden 50 Hyperparameter-Kombinationen getestet, die den
# Hyperparameter-Space so gut wie möglich abdecken.
xgb_res <- 
  xgb_workflow %>% 
  tune_grid(folds,
            grid = 50,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))

# Wir sortieren die Hyperparameter Spezifikationen nach ROC AUC
xgb_res %>% 
  show_best(metric = "roc_auc", n = 10) %>% 
  arrange(desc(mean))

# A tibble: 10 × 9
   trees tree_depth learn_rate .metric .estimator  mean     n std_err .config              
   <int>      <int>      <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
 1  1072          3    0.00597 roc_auc binary     0.890     5  0.0278 Preprocessor1_Model08
 2  1918          3    0.00434 roc_auc binary     0.888     5  0.0285 Preprocessor1_Model06
 3  1255          5    0.0149  roc_auc binary     0.887     5  0.0230 Preprocessor1_Model15
 4   964          2    0.0301  roc_auc binary     0.886     5  0.0188 Preprocessor1_Model05
 5  1158          2    0.0132  roc_auc binary     0.885     5  0.0205 Preprocessor1_Model04
 6  1082         15    0.0228  roc_auc binary     0.884     5  0.0217 Preprocessor1_Model50
 7   562          3    0.00531 roc_auc binary     0.884     5  0.0276 Preprocessor1_Model07
 8  1633          8    0.00677 roc_auc binary     0.884     5  0.0234 Preprocessor1_Model26
 9  1779          8    0.0169  roc_auc binary     0.884     5  0.0217 Preprocessor1_Model27
10   168         13    0.108   roc_auc binary     0.884     5  0.0225 Preprocessor1_Model44

Wir finden so die optimalen Hyperparameterwerte trees = 1072, tree_depth = 3 und learn_rate = 0.00597. Wir können die Tuning-Resultate auch grafisch darstellen:

# Plot Tuning Resultate
autoplot(xgb_res)

Wir speichern das beste Modell ab und verwenden die Cross-Validation Resultate für den Vergleich des XGBoost Modells mit dem Random Forest und dem Classification Tree mittels der ROC Kurven.

# Bestes Modell
xgb_best <- 
  xgb_res %>% 
  select_best(metric = "roc_auc")

# Resultate vorbereiten für ROC Kurve
xgb_auc <- 
  xgb_res %>% 
  collect_predictions(parameters = xgb_best) %>% 
  roc_curve(hd, .pred_yes, event_level = "second") %>% 
  mutate(model = "XGBoost")

# ROC Kurven vergleichen
bind_rows(rf_auc, dt_auc, xgb_auc) %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(linewidth = 1.5, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  scale_color_viridis_d(option = "plasma", end = .6) +
  theme_bw()

Hier ist es tatsächlich so, dass das Random Forest Ensemble am besten funktioniert. Allerdings ist es sehr wahrscheinlich, dass wir beim Tuning des XGBoost Modells noch nicht das Maximum rausgeholt haben. In der Praxis ist typischerweise das XGBoost noch eine Spur besser als der Random Forest.

Wie gewohnt fitten wir am Schluss nun noch das optimale Modell auf dem ganzen Trainingsdatensatz.

# Last Fit Spezifikation (optimale Hyperparameter Werte)
last_xgb_mod <-
  boost_tree(trees = 1072, tree_depth = 3, learn_rate = 0.00597) %>% 
  set_engine("xgboost", event_level = "second", nthread = 8) %>% 
  set_mode("classification")

# Workflow anpassen (optimales Modell)
last_xgb_workflow <- 
  xgb_workflow %>% 
  update_model(last_xgb_mod)

# Last Fit (auf ganzem Trainingsdatensatz)
last_xgb_fit <-
  last_xgb_workflow %>% 
  fit(data = train)

# Modell Output
last_xgb_fit %>% 
  extract_fit_parsnip()

parsnip model object

##### xgb.Booster
raw: 1.2 Mb 
call:
  xgboost::xgb.train(params = list(eta = 0.00597, max_depth = 3, 
    gamma = 0, colsample_bytree = 1, colsample_bynode = 1, min_child_weight = 1, 
    subsample = 1), data = x$data, nrounds = 1072, watchlist = x$watchlist, 
    verbose = 0, nthread = 8, objective = "binary:logistic")
params (as set within xgb.train):
  eta = "0.00597", max_depth = "3", gamma = "0", colsample_bytree = "1", colsample_bynode = "1", min_child_weight = "1", subsample = "1", nthread = "8", objective = "binary:logistic", validate_parameters = "TRUE"
xgb.attributes:
  niter
callbacks:
  cb.evaluation.log()
# of features: 25 
niter: 1072
nfeatures : 25 
evaluation_log:
    iter training_logloss
       1        0.6902110
       2        0.6873069
---                      
    1071        0.1694754
    1072        0.1693849

Auch für dieses Modell können wir uns mithilfe des Packages vip die Wichtigkeit der Variablen anzeigen lassen. Das xgboost Package erstellt intern ein One-Hot-Encoding, weshalb wir hier nun einzelne Ausprägungen von kategorischen Variablen als Variablen sehen. Den mit Abstand grössten Effekt hat die Ausprägung normal der kategorischen Variable thal.

# Variable Importance
last_xgb_fit %>% 
  extract_fit_parsnip() %>% 
  vip(num_features = 20)

Nun vergleichen wir die insgesamt 6 Modelle auf dem Testdatensatz:

# ------------------------------------------------------
# 6. Testset Performance

# Testset-Performance
test_dt_aug <- augment(last_dt_fit, test)
test_rf_aug <- augment(last_rf_fit, test)
test_xgb_aug <- augment(last_xgb_fit, test)

# ROC AUC
test_dt_aug %>% 
  roc_auc(truth = hd, .pred_yes, event_level = "second")
test_rf_aug %>% 
  roc_auc(truth = hd, .pred_yes, event_level = "second")
test_xgb_aug %>% 
  roc_auc(truth = hd, .pred_yes, event_level = "second")

# ACCURACY
test_dt_aug %>% 
  accuracy(truth = hd, .pred_class)
test_rf_aug %>% 
  accuracy(truth = hd, .pred_class)
test_xgb_aug %>% 
  accuracy(truth = hd, .pred_class)

Für die Modelle, die wir nicht mit tidymodels gefittet haben, müssen wir leicht anders vorgehen:

# Performance der Modelle ohne tidymodels
pred_test_bag <- predict(bag_fit, newdata = test, type = "prob")
pred_test_rf <- predict(rf_fit, newdata = test, type = "prob")
pred_test_gbm <- predict(gbm_final_fit, n.trees = gbm_final_fit$n.trees, test, type = "response")

# Resultate in Data Frame Form
test_bag_aug <- data.frame(pred_test_bag) %>% 
  bind_cols(hd = test$hd) %>% 
  rename(.pred_no = no, .pred_yes = yes) %>% 
  mutate(.pred_class = factor(ifelse(.pred_yes >= 0.5, "yes", "no"), levels = c("no", "yes")))
test_rfc_aug <- data.frame(pred_test_rf) %>% 
  bind_cols(hd = test$hd) %>% 
  rename(.pred_no = no, .pred_yes = yes) %>% 
  mutate(.pred_class = factor(ifelse(.pred_yes >= 0.5, "yes", "no"), levels = c("no", "yes")))
test_gbm_aug <- data.frame(pred_test_gbm) %>% 
  bind_cols(hd = test$hd) %>% 
  rename(.pred_yes = pred_test_gbm) %>% 
  mutate(.pred_no = 1 - .pred_yes,
         .pred_class = factor(ifelse(.pred_yes >= 0.5, "yes", "no"), levels = c("no", "yes")))

# ROC AUC
test_bag_aug %>% 
  roc_auc(truth = hd, .pred_yes, event_level = "second")
test_rfc_aug %>% 
  roc_auc(truth = hd, .pred_yes, event_level = "second")
test_gbm_aug %>% 
  roc_auc(truth = hd, .pred_yes, event_level = "second")

# ACCURACY
test_bag_aug %>% 
  accuracy(hd, .pred_class)
test_rfc_aug %>%
  accuracy(hd, .pred_class)
test_gbm_aug %>%  
  accuracy(hd, .pred_class)

Die Resultate sehen wie folgt aus:

Modell	ROC AUC (Testset)	Accuracy (Testset)
Classification Tree (rpart)	83.2%	77.3%
Bagging (randomForest)	90.4%	81.3%
Random Forest (randomForest)	92.9%	86.7%
Random Forest (ranger)	93.3%	90.7%
Gradient Boosting (gbm)	93.0%	89.3%
Gradient Boosting (xgboost)	90.2%	84.0%