Machine Learning mit R und caret: GBM optimieren (Gradient Boosting Machine)

Entscheidungsbaum: Median Eigenheimpreis in Boston

Das Maschinelle Lernen vereinigt Methoden aus unterschiedlichen Fachbereichen. Während Ansätze der klassischen Statistik eher auf Hypothesentests ausgelegt sind, steht beim Data Mining oft die Ableitung von praxisrelevanten Erkenntnissen aus vorhandenen Daten im Vordergrund, und das Machine Learning zielt auf die Anwendung der „trainierten“ Modelle auf zuvor nicht gesehene Daten – sprich Vorhersagen. Bei den jeweils eingesetzten Methoden gibt es jedoch erhebliche Schnittmengen.

Video: GBM optimieren

In diesem Beitrag geht es darum, Preise von Eigenheimen anhand von Charakteristika von Stadtteilen vorherzusagen. Der Datensatz heißt Boston und ist im R-Paket MASS enthalten, das zur Standardinstallation (Base R) gehört. Er kann folgendermaßen aktiviert werden:

library(MASS)       # Paket laden 
data(Boston)        # Daten laden
str(Boston)         # Struktur des Datensatzes betrachten

Die Daten: Preise von Eigenheimen in Boston

Der Datensatz enthält aggregierte Daten aus 506 Stadtteilen (=506 Zeilen / Fälle) von 14 Variablen: dem Median-Wohnungspreis und 13 möglichen Prädiktoren.

Histogramm mit Normalverteilungskurve: Median Wohnungspreise in 1.000 USD in Boston (Darstellung mit dem R-Paket ggplot2)

Die möglichen Einflussfaktoren sind:

  • Kriminalitätsrate pro Kopf (crim)
  • Anteil des Wohnlandes über 25.000 sq.ft (zn)
  • Anteil der Nicht-Handels-Geschäftsfläche (indus)
  • Flusslage am Chas River (chas: 1=am Fluss, 0=nicht am Fluss)
  • Konzentration der Stickstoffoxide (Teile pro 10 Millionen, nox)
  • durchschnittliche Anzahl Zimmer pro Wohnung (rm)
  • Anteil der Eigentumswohnungen, die vor 1940 gebaut wurden (age; der Datensatz wurde 1978 veröffentlicht)
  • gewichtete Entfernung zu fünf Bostoner Beschäftigungszentren (dis)
  • Index der Zugänglichkeit zu Einfallstraßen (rad)
  • vollwertiger Immobilien-Steuersatz pro $10.000 (tax)
  • Schüler-Lehrer-Quotient (ptratio)
  • Anteil Farbiger: (black; Formel: 1000(Bk – 0.63)²)
  • Bevölkerungsanteil mit niedrigem Status (lstat)

Der Einstieg: Entscheidungsbäume

Entscheidungsbäume vereinen Elemente der klassischen Statistik, des Data Mining und des maschinellen Lernens. Ausgehend vom gesamten Datensatz, werden schrittweise Untergruppen gebildet, die sich bei der Zielvariable möglichst deutlich unterscheiden. Dabei können Signifikanzkriterien zum Einsatz kommen (z. B. der Chi-Quadrat-Test) oder andere Maße wie der Komplexitäts-Parameter cp. Im Sinne des Data Mining ist die einfache Interpretierbarkeit: Die Ergebnisse von Entscheidungsbäumen können statistische Laien leicht verstehen. Hier ein rpart-Modell (rpart ist ein Algorithmus und ein R-Paket):

Entscheidungsbaum: Median Eigenheimpreis in Boston
Entscheidungsbaum: Median Eigenheimpreis in 1.000 US-Dollar in Boston (Klicken für größere Ansicht)

Im Sinne der Lesbarkeit wurde das Modell vereinfacht – weitere Verzweigungen unterhalb der zweiten Ebene wären möglich.

  • Der durchschnittliche Eigenheimpreis des gesamten Datensatzes beträgt rund 23.000 USD (Knoten 1 oben).
  • Die deutlichste Differenzierung ergibt sich nach Anzahl der Räume (Durchschnitt pro Stadtteil, rm): Unter 6,9 (20.000 USD, Verzweigung nach links) vs. mindestens 6,9 (37.000 USD, Verzweigung nach rechts).
  • Im rechten Ast wird erneut nach der Anzahl der Räume aufgeteilt. Der höchste Median-Eigenheimpreis in diesem Modell liegt im Knoten 7, unten rechts: 45.000 USD bei durchschnittlich 7,4 Räumen im Stadtteil.
  • Der geringste Median-Eigenheimpreis betrifft Stadtteile mit durchschnittlich weniger als 6,9 Räumen pro Einheit mit einem Bevölkerungsanteil mit niedrigem Status (lstat) von mindestens 14% (links unten, Knoten 3, 15.000 USD).

Der Machine Learning-Anteil besteht darin, dass es der Algorithmus dem Anwender erspart, selbst nach Verzweigungen zu suchen, d. h. verschiedene Merkmale manuell auf Unterschiede bei der Zielvariable zu testen. Desweiteren übernimmt der Algorithmus auch die Entscheidung, Trennwerte bei kontinuierlichen Merkmalen zu finden – hier: bei welchem Wert werden Untergruppen nach Anzahl der Räume oder Bevölkerungsanteil mit niedrigem Status gebildet. Dies ist eine praktische Alternative zu den früher üblichen Tabellenbänden, in denen man lange nach der „Nadel im Heuhaufen“ suchen konnte – oder, bei ungünstiger Gruppenbildung bei kontinuierlichen Merkmalen (typisch z. B. Altersgruppen), relevante Unterschiede verpasste.

Die Herausforderung: GBM vs. Random Forest

Nehmen wir an, wir nehmen an einem Machine-Learning-Wettbewerb teil. Die Aufgabe: Eine möglichst hohe Prognose-Genauigkeit mit einem Modell zu erreichen, das Eigenheimpreise auf Basis der oben genannten Merkmale von Stadtteilen vorhersagt. Nehmen wir weiter an, ein Kollege hat schon vorgelegt: Er hat ein Random-Forest-Modell aufgestellt. Random Forests basieren auf einer Menge an Entscheidungsbäumen – die finale Prognose wird aus den vielen Bäumen gemittelt. Der Trick besteht darin, unterschiedliche Bäume aufzustellen, die gegenseitig ihre Schwächen ausgleichen können. Neben den unterschiedlichen Zufallsstichproben wird das zusätzlich dadurch erreicht, dass bei jeder Verzweigung nicht alle Prädiktoren zur Auswahl stehen, sondern nur eine zufällig gewählte Teilmenge. (Klingt absurd? Das klang auch für führende Profis anfangs so, doch die Technik hat sich in vielen Anwendungen bewährt!)

Wir versuchen, mit einem GBM-Modell (Gradient Boosting Machine; das R-Paket gbm spricht von Generalized Boosted Regression Modeling, da die Grundidee auf unterschiedliche Algorithmen verallgemeinert wird) zu kontern. Ausgangspunkt sind auch hier Entscheidungsbäume. Allerdings wird hier jeder Baum auf die Residuen (Vorhersagefehler) des vorherigen Baums angesetzt, um so schrittweise die Abweichungen zwischen Modell und beobachteten Zieldaten zu minimieren. Dies geschieht, indem die schlecht modellierten Fälle stärker gewichtet werden. Im Gegensatz zum Random Forest-Modell sind die Bäume damit nicht von einander unabhängig. Das Motto des Boostings: Mache einen „schwachen Lerner“ [Entscheidungsbaum] schrittweise zu einem starken Lerner.

Kreuzvalidierung

Wie gut sind unsere Modelle geeignet, Eigenheimpreise bei neuen Daten anhand der genannten Merkmale vorherzusagen? Da uns keine neuen Daten zur Verfügung stehen, wenden wir einen Trick aus dem Machine Learning an: Wir simulieren Trainings- und Testdaten, indem wir wiederholt Zufallsstichproben aus den Daten ziehen. Ein Teil der Daten wird verwendet, um das Modell zu schätzen (Trainingsdaten), während ein anderer Teil zurückgehalten wird, um das Modell auf seine Vorhersagegenauigkeit zu testen (Testdaten). Da das Testergebnis je nach Aufteilung besser oder schlechter ausfallen kann, wird dieser Vorgang öfters wiederholt.

5-fache Kreuzvalidierung (Schema)
5-fache Kreuzvalidierung (Schema)

Wir verwenden hier 10-fache Kreuzvalidierung (im Bild ist der Einfachheit halber nur 5-fache Kreuzvalidierung dargestellt): Die Daten werden in 10 gleich große Blöcke aufgeteilt, jeweils ein Teil bildet den Testdatensatz, die anderen 9 Teile die Trainingsdaten. Dieser Vorgang wird 10 mal wiederholt, bis jedes Zehntel ein mal als Testdatensatz fungierte. Schließlich machen wir uns die moderne Rechenpower zunutze, indem wir diesen Vorgang insgesamt 10 mal wiederholen, mit unterschiedlichen Aufteilungen in 10 Teile. So erhalten wir recht stabile Schätzungen der Vorhersage-Genauigkeit.

Kreuzvalidierung leicht gemacht: Das caret-Paket

Wer möchte, kann die Kreuzvalidierung manuell programmieren. Einfacher ist es, das caret-Paket von Max Kuhn zu nutzen, dem Autor von Applied Predictive Modeling. Im folgenden vergleichen wir einen einzelnen Entscheidungsbaum mit dem Random Forest-Modell des Kollegen.

caret bietet eine Schnittstelle zu einer Vielzahl von Machine-Learning-Algorithmen. Wir greifen auf die train-Funktion zu und müssen uns nicht mit der unterschiedlichen Syntax von Machine-Learning-Methoden befassen.

trainctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 10)
set.seed(2018)
rpart_tree <- train(medv ~ ., data = Boston, method = "rpart", trControl = trainctrl)
rf_tree <- train(medv ~ ., data = Boston, method = "rf",
 trControl = trainctrl)
  • Zunächst wird die Kreuzvalidierungsmethode definiert: hier 10-fache Kreuzvalidierung mit 10 Wiederholungen.
  • Der Zufallsgenerator wird auf einen festen Wert gesetzt (set.seed), um die Ergebnisse reproduzierbar zu machen – so werden bei wiederholter Ausführung die gleichen Teildatensätze für die Kreuzvalidierung gebildet.
  • Sodann wird das Modell mit einem einzelnen Entscheidungsbaum aufgestellt: rpart_tree. Im Vergleich zu einem Entscheidungsbaum ohne Kreuzvalidierung wird hier bereits Overfitting (Überanpassung) minimiert.
  • Das Random-Forest-Modell (rf_tree) wird über die gleiche Syntax erstellt – mit dem einen Unterschied, dass wir „rf“ als Methode spezifizieren.
    (Wir nehmen an, dass wir dieses Modell vom Kollegen erhalten haben – mit einem triumphierenden Blick und der Gewissheit, Random Forest sei hier nicht zu schlagen. Wir nehmen die Herausforderung dankend an.)
  • Die Schreibweise „medv ~ .“ ist die sog. Formelschreibweise, die in R häufig verwendet wird. Vor der Tilde (~) steht die abhängige Variable. Der Punkt ist die kürzeste Art, R mitzuteilen, dass alle anderen Variablen im Datensatz zur Modellierung verwendet werden sollen.

Nun wollen wir die Kreuzvalidierung durchführen und die Ergebnisse vergleichen:

resamps <- resamples(list(rpart_tree = rpart_tree, randomForest = rf_tree))
summary(resamps)

Dazu greifen wir auf carets resamples-Funktion zurück und übergeben ihr einfach eine Liste mit den zu vergleichenden Modellen. Mit der summary-Funktion erhalten wir eine numerische Zusammenfassung. Als Gütekriterium verwenden wir RMSE, den root mean squared error (Wurzel des quadrierten Vorhersagefehlers).  Durch die Quadrierung ist sichergestellt, dass sich positive und negative Vorhersagefehler nicht aufheben. RMSE hat den Charme, in Einheiten der abhängigen Variable interpretierbar zu sein:

  • Der einzelne Entscheidungsbaum kommt auf einen durchschnittlichen RMSE (Durchschnitt über die 10 Wiederholungen von 10 Kreuzvalidierungen, also 100 Evaluierungsdurchgängen) von 5,68. Das entspricht einem durchschnittlichen Fehler von 5.680 US-Dollar (eine Einheit entspricht 1.000 USD).
  • Das Random-Forest-Modell ist mit 3.09 deutlich genauer, d. h. es weicht mit den modellierten Werten durchschnittlich nur um 3.090 US-Dollar von den tatsächlich beobachteten Werten ab.

GBM vs. Random Forest

Wie schlägt sich unser Gradient Boosting Machine-Modell im Vergleich?

set.seed(2018)
gbm_tree_auto <- train(medv ~ ., data = Boston, method = "gbm", distribution = "gaussian", trControl = trainctrl, verbose = FALSE)
gbm_tree_auto

Wir verwenden die gleiche Kreuzvalidierungsmethode wie oben, die im Objekt trainctrl gespeichert ist (10fach mit 10 Wiederholungen). Es folgt wieder ein Zugriff auf die train-Funktion, diesmal mit der Methode „gbm“. verbose = FALSE unterdrückt eine sehr umfangreiche Ausgabe von Ergebnissen. Wir geben die Gauss-Verteilung an, da es sich um ein metrisches Merkmal handelt.

caret stellt nicht nur ein Gradient Boosting-Machine-Modell auf, sondern optimiert automatisch anhand der Kreuzvalidierungsergebnisse einige Modellparameter.

Enttäuschend: Das gbm-Modell schlägt zwar erwartungsgemäß den einzelnen Entscheidungsbaum sehr deutlich, bleibt mit einem RMSE von 3,20 jedoch hinter dem Random-Forest-Modell zurück.

GBM und caret: Modell-Optimierung

So schnell wollen wir uns jedoch dem triumphierenden Kollegen nicht geschlagen geben. Was können wir tun? Lesen wir die Ausgabe zu unserem gbm_tree_auto-Modell von oben genauer, so sehen wir, dass caret vier Modellparameter festgelegt hat. Diese Parameter können wir auch über modelLookup(„gbm“) abrufen:

  • n.trees: Die Anzahl der Entscheidungsbäume. Caret prüfte 50, 100 und 150 und entschied sich für 150.
  • interaction.depth: Die „Tiefe“ der Bäume, d. h. die Anzahl der Ebenen. (Der Beispiel-Baum ganz oben enthält zwei Verzweigungsebenen.) Caret prüfte 1, 2 und 3 und entschied sich für 3.
  • shrinkage: Die „Lerngeschwindigkeit“ – betrifft die Gewichtung schlecht vorhergesagter Fälle für den nachfolgenden Baum. Caret fixierte diesen Parameter auf 0,1.
  • n.minobsinnode: Die Mindestzahl an Beobachtungen (Fällen) pro Knoten, d. h. die Mindestgröße der Knoten. Caret fixierte diesen Wert auf 10.

Die Kreuzvalidierung ist sehr rechenaufwändig. Caret trifft daher intelligente Vorannahmen über die Modellparameter. Es ist jedoch nicht garantiert, dass diese zum bestmöglichen Modell führen. Wir können versuchen, andere Werte für die Modellparameter vorzugeben und hoffen, dadurch das Modell zu verbessern. Praktisch: Wir müssen nicht einzeln jede Kombination prüfen, die wir für sinnvoll halten, sondern können ein ganzes Raster an Kombinationen vorgeben, die caret systematisch abarbeitet und per Kreuzvalidierung bewertet. Zum Beispiel so:

myGrid <- expand.grid(n.trees = c(150, 175, 200, 225),
 interaction.depth = c(5, 6, 7, 8, 9),
 shrinkage = c(0.075, 0.1, 0.125, 0.15, 0.2),
 n.minobsinnode = c(7, 10, 12, 15))

set.seed(2018)
gbm_tree_tune <- train(medv ~ ., data = Boston, method = "gbm", distribution = "gaussian",
 trControl = treectrl, verbose = FALSE,
 tuneGrid = myGrid)
gbm_tree_tune

  • Wir legen ein Raster (grid) fest. Die Base-R-Funktion expand.grid bildet eine Matrix mit allen Kombinationen der vier Parameter.
  • Dieses Raster übergeben wir mittels des tuneGrid-Parameters der train-Funktion.
  • Nun brauchen wir etwas Geduld … bis das Ergebnis kommt:
    • n.trees = 200 (statt zuvor 150)
    • interaction.depth = 9 (statt zuvor 3)
    • shrinkage = 0.1 (wie oben)
    • n.minobsinnode = 7 (statt zuvor 10)

Mit diesen Parametern gehen wir nun erneut in den Wettbewerb – wie fällt der Vergleich mit dem Random-Forest-Modell jetzt aus?

myGrid <- gbm_tree_tune$bestTune

set.seed(2018)
gbm_tree <- train(medv ~ ., data = Boston, method = "gbm",
 trControl = trainctrl,
 tuneGrid = myGrid, verbose = FALSE)
gbm_tree

resamps <- resamples(list(rpart_tree = rpart_tree, randomForest = rf_tree, gbm_tree_auto = gbm_tree_auto,gbm_tree_user = gbm_tree))
summary(resamps)
  • Oben haben wir ein tuning ausgeführt, das Objekt gbm_tree_tune steht uns nun zur Verfügung. Auf die Parameter des optimierten Modells können wir mittels gbm_tree_tune$bestTune zugreifen und diese als Raster festlegen. Im Folgenden wird dann nur dieses eine Modell für den Modellvergleich verwendet. In der Ausgabe (gbm_tree) finden wir die optimierten Modellparameter wieder.
  • Wie oben vergleichen wir mittels resamples die Modelle.
  • Ergebnis: Diesmal gewinnen wir – mit einem durchschnittlichen RMSE von 3.088 ist unser GBM-Modell minimal besser als das Random-Forest-Modell (RMSE = 3.091).

Grafisch:

dotplot(resamps, metric = "RMSE", main = "Modell-Vergleich")
caret: Modellvergleich dotplot
Modellvergleich mit caret: dotplot. Random Forest vs. GBM (gbm_tree_auto), GBM mit optimierten Parametern (gbm_tree_user) und einzelnem Entscheidungsbaum (rpart_tree)

caret ordnet die Modelle in der Reihenfolge der Kennzahl (hier: RMSE) an. Unser benutzerdefiniertes GBM-Modell schneidet jetzt tatsächlich minimal besser ab als das Random-Forest-Modell. An dritter Stelle findet sich das GBM-Modell mit den Parametern ein, die caret in der Voreinstellung gefunden hat. Der einzelne Entscheidungsbaum ist wie erwartet weit abgeschlagen.

Zusammenfassung: Modelloptimierung mit caret

Der Beitrag zeigte, wie man mit Hilfe des Pakets caret und einer einheitlichen Schnittstelle (train-Funktion) verschiedene Modelle (einzelner Entscheidungsbaum mit rpart, Random Forest, GBM) aufstellen und per Kreuzvalidierung vergleichen kann. Mittels eigener Vorgaben zu den Modellparametern gelang es, das ursprünglich schwächere GBM-Modell so zu optimieren, dass es das Random-Forest-Modell knapp übertraf.

Manche Autoren empfehlen, nicht das „beste“ Modell zu nehmen, sondern ein einfacheres, das nicht weiter als eine Standardabweichung vom besten Modell abweicht („one standard error rule“ bzw. „1 SE rule“). Dies ist eine zusätzliche Maßnahme, um Überanpassung (overfitting) zu vermeiden.

Nicht probiert haben wir, das Random-Forest-Modell seinerseits gegenüber der caret-Voreinstellung zu optimieren. Hier steht nur ein Parameter zur Auswahl: Die Anzahl der Prädiktoren, die für jede Verzweigung per Zufallsauswahl herangezogen werden (mtry).

In einem Machine Learning-Workshop steht mehr Zeit zur Verfügung, um diese und weitere Techniken kennen zu lernen und zu üben.

Freue mich über Kommentare!