6  tidymodels

6.1 Lernsteuerung

6.1.1 Lernziele

  • Sie sind in der Lage, Regressionsmodelle mit dem tidymodels-Ansatz zu spezifizieren.
  • Sie können Begriffe des statistischen Lernens in das Vokabular von tidymodels übersetzen.

6.2 Vorbereitung

  • Lesen Sie TMWR, Kapitel 1
  • Lesen Sie übrige Literatur zu diesem Thema: TMWR, Kap. 1, 5, 6, 7, 8, 9

6.2.1 Benötigte R-Pakete

tidymodels ist ein Metapaket: Ein (R-)Paket, das mehrere andere Paket startet und uns damit das Leben einfacher macht, analog zu tidyverse. Eine Liste der R-Pakete, die durch tidymodels gestartet werden, findet sich hier. Probieren Sie auch mal ?tidymodels.

Eine Liste aller Pakete, die in Tidymodels benutzt werden, die dependencies, kann man sich so ausgeben lassen:

pkg_deps(x = "tidymodels", recursive = FALSE)

6.3 Daten

Dieser Abschnitt bezieht sich auf Kapitel 4 in Silge und Kuhn (2022).

Wir benutzen den Datensatz zu Immobilienpreise aus dem Ames County in Iowa, USA, gelegen im Zentrum des Landes.

data(ames)  # Daten wurden über tidymodels mit geladen
ames <- 
  ames %>% 
  mutate(Sale_Price = log10(Sale_Price))

Hier wurde die AV log-transformiert. Das hat zwei (wichtige) Effekte:

  1. Die Verteilung ist symmetrischer, näher an der Normalverteilung. Damit gibt es mehr Daten im Hauptbereich des Ranges von Sale_Price, was die Vorhersagen stabiler machen dürfte.
  2. Logarithmiert man die Y-Variable, so kommt dies einem multiplikativen Modell gleich, s. auch hier.

6.4 Train- vs Test-Datensatz aufteilen

Dieser Abschnitt bezieht sich auf Kapitel 5 in Silge und Kuhn (2022).

Hinweis

Das Aufteilen in Train- und Test-Datensatz ist einer der wesentlichen Grundsätze im maschinellen Lernen. Das Ziel ist, Overfitting abzuwenden. Im Train-Datensatz werden alle Modelle berechnet. Der Test-Datensatz wird nur einmal verwendet, und zwar zur Überprüfung der Modellgüte.

Eine Faustregel ist es, 70-80% der Daten in das Train-Sample und die übrigen 20-30% in das Test-Sample zu stecken, s. Abbildung 6.1

pie title Train-Test-Aufteilung
    "Train" : 80
    "Test" : 19
    "The Unkown God": 1

Abbildung 6.1: 80-20-Aufteilung der Daten in Train- bzw. Test-Sample

Praktisch funktioniert das in Silge und Kuhn (2022) wie folgt.

Wir laden die Daten und erstellen einen Index, der jeder Beobachtung die Zuteilung zu Train- bzw. zum Test-Datensatz zuweist.

Das kann, mit tidymodels so aussehen:

ames_split <- initial_split(ames, prop = 0.80, strata = Sale_Price)

initial_split() speichert für spätere komfortable Verwendung auch die Daten. Aber eben auch der Index, der bestimmt, welche Beobachtung im Train-Set landet:

ames_split$in_id %>% head(n = 10)
##  [1] 28 30 31 32 33 35 79 83 84 87
length(ames_split$in_id)
## [1] 2342

Praktisch ist auch, dass die AV-Verteilung in beiden Datensätzen ähnlich gehalten wird (Stratifizierung), das besorgt das Argument strata.

Die eigentlich Aufteilung in die zwei Datensätze geht dann so:

ames_train <- training(ames_split)
ames_test  <-  testing(ames_split)

6.5 Grundlagen der Modellierung mit tidymodels

Dieser Abschnitt bezieht sich auf Kapitel 6 in Silge und Kuhn (2022).

tidymodels ist eine Sammlung mehrerer, zusammengehöriger Pakete, eben zum Thema statistische Modellieren.

Das kann man analog zur Sammlung tidyverse verstehen, zu der z.B. das R-Paket dplyr gehört.

Das R-Paket innerhalb von tidymodels, das zum “Fitten” von Modellen zuständig ist, heißt parsnip.

Eine Liste der verfügbaren Modelltypen, Modellimplementierungen und Modellparameter, die in Parsnip aktuell unterstützt werden, findet sich hier.

6.5.1 Modelle spezifizieren

Ein (statistisches) Modell wird in Tidymodels mit drei Elementen spezifiziert, vgl. Abbildung 6.2.

flowchart LR
   
  

  subgraph Modus
  r2[regresssion]
  classification
  end
  
  subgraph Implementierung
  lm
  stan_glm
  div2[...]
  end
  
  subgraph Algorithmus
  R[Regression]
  NN[Neuronale Netze]
  div[...]
  end

Abbildung 6.2: Definition eines Models in tidymodels

Die Definition eines Modells in tidymodels folgt diesen Ideen:

  1. Das Modell sollte unabhängig von den Daten spezifiziert sein
  2. Das Modell sollte unabhängig von den Variablen (AV, UVs) spezifiziert sein
  3. Das Modell sollte unabhängig von etwaiger Vorverarbeitung (z.B. z-Transformation) spezifiziert sein

Da bei einer linearen Regression nur der Modus “Regression” möglich ist, muss der Modus in diesem Fall nicht angegeben werden. Tidymodels erkennt das automatisch.

lm_model <-   
  linear_reg() %>%   # Algorithmus, Modelltyp
  set_engine("lm")  # Implementierung
  # Modus hier nicht nötig, da lineare Modelle immer numerisch klassifizieren

6.5.2 Modelle berechnen

Nach Rhys (2020) ist ein Modell sogar erst ein Modell, wenn die Koeffizienten berechnet sind. Tidymodels kennt diese Unterscheidung nicht. Stattdessen spricht man in Tidymodels von einem “gefitteten” Modell, sobald es berechnet ist. Ähnlich fancy könnte man von einem “instantiierten” Modell sprechen.

Für das Beispiel der einfachen linearen Regression heißt das, das Modell ist gefittet, sobald die Steigung und der Achsenabschnitt (sowie die Residualstreuung) berechnet sind.

lm_form_fit <- 
  lm_model %>% 
  fit(Sale_Price ~ Longitude + Latitude, data = ames_train)

6.5.3 Vorhersagen

Im maschinellen Lernen ist man primär an den Vorhersagen interessiert, häufig nur an Punktschätzungen. Schauen wir uns also zunächst diese an.

Vorhersagen bekommt man recht einfach mit der predict() Methode von tidymodels1:

predict(lm_form_fit, new_data = ames_test) %>% 
  head()

Die Syntax zum Vorhersagen lautet also: predict(modell, daten_zum_vorhersagen).

6.5.4 Vorhersagen im Train-Datensatz

Vorhersagen im Train-Datensatz machen kaum Sinn, da sie nicht gegen Overfitting geschützt sind und daher deutlich zu optimistisch sein können.

Bei einer linearen Regression ist diese Gefahr nicht so hoch, aber bei anderen, flexibleren Modellen, ist diese Gefahr absurd groß.

6.5.5 Modellkoeffizienten im Train-Datensatz

Gibt man den Namen des Modellobjekts ein, so wird ein Überblick an relevanten Modellergebnissen am Bildschirm gedruckt:

lm_form_fit
## parsnip model object
## 
## 
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
## 
## Coefficients:
## (Intercept)    Longitude     Latitude  
##    -314.955       -2.134        2.863

Innerhalb des Ergebnisobjekts findet sich eine Liste namens fit, in der die Koeffizienten (der “Fit”) abgelegt sind:

lm_form_fit %>% pluck("fit")
## 
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
## 
## Coefficients:
## (Intercept)    Longitude     Latitude  
##    -314.955       -2.134        2.863

Zum Herausholen dieser Infos kann man auch alternativ die Funktion extract_fit_engine() verwenden:

lm_fit <-
  lm_form_fit %>% 
  extract_fit_engine()

lm_fit
## 
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
## 
## Coefficients:
## (Intercept)    Longitude     Latitude  
##    -314.955       -2.134        2.863
Hinweis

Möchten Sie wissen, was sich in lm_form_fit alles verbirgt, bietet sich die Funktion str an. Alternativ können Sie in RStudio unter Environment das Objekt “aufklappen”.

Das extrahierte Objekt ist, in diesem Fall, das typische lm() Objekt. Entsprechend kann man daruaf coef() oder summary() anwenden.

coef(lm_fit)
## (Intercept)   Longitude    Latitude 
## -314.954598   -2.133832    2.863270
summary(lm_fit)
## 
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02404 -0.09525 -0.01574  0.09584  0.54193 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -314.9546    14.3978  -21.88   <2e-16 ***
## Longitude     -2.1338     0.1274  -16.75   <2e-16 ***
## Latitude       2.8633     0.1804   15.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1588 on 2339 degrees of freedom
## Multiple R-squared:  0.1794, Adjusted R-squared:  0.1787 
## F-statistic: 255.6 on 2 and 2339 DF,  p-value: < 2.2e-16

Schicker sind die Pendant-Befehle aus broom, die jeweils einen Tibble zuückliefern:

library(broom)
tidy(lm_fit) # Koeffizienten
glance(lm_fit) # Modellgüte

Eine weitere Alternative sind die Befehle zur Modell-Performance von easystats´^[Paketperformance`]:

library(easystats)
parameters(lm_form_fit)
r2(lm_form_fit)
## # R2 for Linear Regression
##        R2: 0.179
##   adj. R2: 0.179
mae(lm_form_fit)
## [1] 0.1211582

6.5.6 Parsnip RStudio add-in

Mit dem Add-in von Parsnip kann man sich eine Modellspezifikation per Klick ausgeben lassen. Nett!

parsnip_addin()

6.6 Workflows

Dieser Abschnitt bezieht sich auf Kapitel 7 in Silge und Kuhn (2022).

6.6.1 Konzept des Workflows in Tidymodels

Definition eines Models in tidymodels

6.6.2 Einfaches Beispiel

Wir initialisieren einen Workflow, verzichten auf Vorverarbeitung und fügen ein Modell hinzu:

lm_workflow <- 
  workflow() %>%  # init
  add_model(lm_model) %>%   # Modell hinzufügen
  add_formula(Sale_Price ~ Longitude + Latitude)  # Modellformel hinzufügen

Werfen wir einen Blick in das Workflow-Objekt:

lm_workflow
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude + Latitude
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear Regression Model Specification (regression)
## 
## Computational engine: lm

Wie man sieht, gehört die Modellformel (y ~ x) zur Vorverarbeitung aus Sicht von Tidymodels.

Was war nochmal im Objekt lm_model enthalten?

lm_model
## Linear Regression Model Specification (regression)
## 
## Computational engine: lm

Jetzt können wir das Modell berechnen (fitten):

lm_fit <- 
  lm_workflow %>%
  fit(ames_train)

Natürlich kann man synonym auch schreiben:

lm_fit <- fit(lm_wflow, ames_train)

Schauen wir uns das Ergebnis an:

lm_fit
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude + Latitude
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:
## stats::lm(formula = ..y ~ ., data = data)
## 
## Coefficients:
## (Intercept)    Longitude     Latitude  
##    -314.955       -2.134        2.863

6.6.3 Vorhersage mit einem Workflow

Die Vorhersage mit einem Tidymodels-Workflow ist einerseits komfortabel, da man einfach sagen kann:

“Nimm die richtigen Koeffizienten des Modells aus dem Train-Set und wende sie auf das Test-Sample an. Berechne mir die Vorhersagen und die Modellgüte.”

So sieht das aus:

final_lm_res <- last_fit(lm_workflow, ames_split)
final_lm_res

Also, last_fit kümmert sich um Folgendes:

  1. Berechne Modell im (kompletten) Train-Sample
  2. Sage Daten im Test-Sample vorher
  3. Berechne Modellgüte im Test-Sample

Es wird ein recht komplexes Objekt zurückgeliefert, das man erst mal durchschauen muss.

Wie man sieht, gibt es mehrere Listenspalten in final_lm_res. Besonders interessant erscheinen natürlich die Listenspalten .metrics und .predictions.

Schauen wir uns die Vorhersagen an. Diese finden sich im resultierenden Objekt von last_fit, zusammen mit anderen Informationen wie MOdellgüte. Die .predictions sind selber ein Tibble, wo in der ersten Spalte die Vorhersagen stehen.

lm_preds <- final_lm_res %>% pluck(".predictions", 1)

Es gibt auch eine Funktion, die obige Zeile vereinfacht (also synonym ist):

lm_preds <- collect_predictions(final_lm_res)
lm_preds %>% slice_head(n = 5)

6.6.4 Modellgüte

Dieser Abschnitt bezieht sich auf Kapitel 9 in Silge und Kuhn (2022).

Die Vorhersagen bilden die Basis für die Modellgüte (“Metriken”), die schon fertig berechnet im Objekt final_lm_res liegen und mit collect_metrics herausgenommen werden können:

lm_metrics <- collect_metrics(final_lm_res)

Alternativ kommt man mit pluck(final_lm_res, ".metrics") an die gleichen Informationen.

.metric .estimator .estimate .config
rmse standard 1.70 × 10−1 Preprocessor1_Model1
rsq standard 1.48 × 10−1 Preprocessor1_Model1

Man kann auch angeben, welche Metriken der Modellgüte man bekommen möchte:

ames_metrics <- metric_set(rmse, rsq)

ames_metrics(data = lm_preds, 
             truth = Sale_Price, 
             estimate = .pred)

6.6.5 Vorhersage von Hand

Man kann sich die Metriken auch von Hand ausgeben lassen, wenn man direktere Kontrolle haben möchte als mit last_fit und collect_metrics.

ames_test_small <- ames_test %>% slice(1:5)
predict(lm_form_fit, new_data = ames_test_small)

Jetzt binden wir die Spalten zusammen, also die “Wahrheit” (\(y\), die beobachteten, tatsächlichen Y-Werte) und die Vorhersagen (\(\hat{y}\)):

ames_test_small2 <- 
  ames_test_small %>% 
  select(Sale_Price) %>% 
  bind_cols(predict(lm_form_fit, ames_test_small)) %>% 
  # Add 95% prediction intervals to the results:
  bind_cols(predict(lm_form_fit, ames_test_small, type = "pred_int")) 
rsq(ames_test_small2, 
   truth = Sale_Price,
   estimate = .pred
   )

Andere Koeffizienten der Modellgüte können mit rmse oder mae2 abgerufen werden.

6.7 Rezepte zur Vorverarbeitung

Dieser Abschnitt bezieht sich auf Kapitel 8 in Silge und Kuhn (2022).

6.7.1 Was ist Rezept und wozu ist es gut?

So könnte ein typischer Aufruf von lm() aussehen:

lm(Sale_Price ~ Neighborhood + log10(Gr_Liv_Area) + Year_Built + Bldg_Type, 
   data = ames)

Neben dem Fitten des Modells besorgt die Formel-Schreibweise noch einige zusätzliche nützliche Vorarbeitung:

  1. Definition von AV und AV
  2. Log-Transformation von Gr_Liv_Area
  3. Transformation der nominalen Variablen in Dummy-Variablen

Das ist schön und nützlich, hat aber auch Nachteile:

  1. Das Modell wird nicht nur spezifiziert, sondern auch gleich berechnet. Das ist unpraktisch, weil man die Modellformel vielleicht in anderen Modell wiederverwenden möchte. Außerdem kann das Berechnen lange dauern.
  2. Die Schritte sind ineinander vermengt, so dass man nicht einfach und übersichtlich die einzelnen Schritte bearbeiten kann.

Praktischer wäre also, die Schritte der Vorverarbeitung zu ent-flechten. Das geht mit einem “Rezept” aus Tidymodels:

simple_ames <- 
  recipe(Sale_Price ~ Neighborhood + Gr_Liv_Area + Year_Built + Bldg_Type,
         data = ames_train) %>%
  step_log(Gr_Liv_Area, base = 10) %>% 
  step_dummy(all_nominal_predictors())
simple_ames
Hinweis

Ein Rezept berechnet kein Modell. Es macht nichts außer die Vorverarbeitung des Modells zu spezifizieren (inklusive der Modellformel).

6.7.2 Workflows mit Rezepten

Jetzt definieren wir den Workflow nicht nur mit einer Modellformel, sondern mit einem Rezept:

lm_workflow <-
  workflow() %>% 
  add_model(lm_model) %>% 
  add_recipe(simple_ames)

Sonst hat sich nichts geändert.

Wie vorher, können wir jetzt das Modell berechnen und uns im Test-Set die Vorhersagen berechnen lassen:

final_lm_res <- last_fit(lm_workflow, ames_split)
final_lm_res

Hier ist die Modellgüte:

lm_metrics <- collect_metrics(final_lm_res)
lm_metrics

6.7.3 Spaltenrollen

Eine praktische Funktion ist es, bestimmte Spalten nicht als Prädiktor, sondern als ID-Variable zu nutzen. Das kann man in Tidymodels komfortabel wie folgt angeben:

ames_recipe <-
  simple_ames %>% 
  update_role(Neighborhood, new_role = "id")

ames_recipe

6.7.4 Preppen und Backen

Ein Rezept ist erstmal nur, ja, ein Rezept: Eine Beschreibung von Schritten und Zutaten. Es ist noch nichts gebacken. Um aus einen Rezept einen “Kuchen” - den transformierten Datensatz - zu bekommen, sind zwei Schritte nötig:

  1. Vorbereiten (to prep): Die Parameter des Rezeptschritte berechnen. So muss der Schritt step_center(var) den Mittelwert von var wissen, sonst kann der Schritt nicht durchgeführt werden.
  2. Backen ist das Rezept vorbereitet, kann der Datensatz damit gebacken werden.

Praktischerweise erledigt Tidymodels das alles automatisch für uns, wir haben da nichts zu tun.

Allerdings ist es manchmal praktisch, den durch das Rezept “gebackenen” (transformierten) Datensatz zu sehen, daher sollte man wissen, wie man das “preppen” und “backen” von Hand erledigt.

  1. Preppen:
ames_recipe_prepped <-
  ames_recipe %>% 
  prep()

ames_recipe_prepped
  1. Backen:
ames_train_baked <- 
  ames_recipe_prepped %>% bake(new_data = NULL) 

ames_train_baked %>% 
  head()

6.7.5 Fazit

Mehr zu Rezepten findet sich hier. Ein Überblick zu allen Schritten der Vorverarbeitung findet sich hier.

6.8 Aufgaben

  1. tidymodels-ames-01
  2. tidymodels-ames-02
  3. tidymodels-ames-03
  4. tidymodels-ames-04
  5. bike01

6.9 Fallstudien

  1. im Gegensatz zum predict() von lm mit Unterstrich bei new_data, also nicht newdata.↩︎

  2. Achtung: Die Funktion mae gibt es sowohl in tidymodels auch in easystats, hier kann es zu Konflikten kommen.↩︎