Abstimmung von Random Forest auf Zeitreihendaten

„Schon wieder Steuern und Random Forest?“ Mein Kollege Thomas hier bei statworx zog die Augenbrauen hoch, als ich ihm von diesem Beitrag erzählte. „Ja!“, antwortete ich, „aber nicht, weil ich Steuern so sehr liebe (wer tut das schon?). Diesmal geht es ums Tuning!“

Spulen wir kurz zurück: Im vorherigen Beitrag haben wir uns angesehen, wie wir ökonometrische Techniken wie Differenzierung mit Machine-Learning-Algorithmen wie Random Forest kombinieren können, um eine Zeitreihe mit hoher Genauigkeit vorherzusagen. Falls du ihn verpasst hast, empfehle ich dir, ihn dir hier anzusehen. Die Daten sind inzwischen auch als CSV-Datei auf unserem statworx-GitHub verfügbar.

Da wir im letzten Beitrag schon einiges abgedeckt haben, blieb kaum Platz für andere Themen. Genau hier setzt dieser Beitrag an. Heute schauen wir uns an, wie wir die Hyperparameter eines Random Forests beim Umgang mit Zeitreihendaten tunen können.

Interessiert? Alles klar, dann legen wir los! Wenn du den letzten Beitrag gelesen hast, kannst du Abschnitt 1 überspringen und direkt mit Abschnitt 2 weitermachen.

Das Setup

# load packages
suppressPackageStartupMessages(require(tidyverse))
suppressPackageStartupMessages(require(tsibble))
suppressPackageStartupMessages(require(randomForest))
suppressPackageStartupMessages(require(forecast))

# specify the path to the csv file (your path here)
file <- "/Users/manueltilgner/Library/Mobile Documents/com~apple~CloudDocs/Artificial Intelligence/10_Blog/RF/tax.csv"

# read in the csv file
tax_tbl <- readr::read_delim(
  file = file,
  delim = ";",
  col_names = c("Year", "Type", month.abb),
  skip = 1,
  col_types = "iciiiiiiiiiiii",
  na = c("...")
) %>% 
  select(-Type) %>% 
  gather(Date, Value, -Year) %>% 
  unite("Date", c(Date, Year), sep = " ") %>% 
  mutate(
    Date = Date %>% 
      lubridate::parse_date_time("m y") %>% 
      yearmonth()
  ) %>% 
  drop_na() %>% 
  as_tsibble(index = "Date") %>% 
  filter(Date <= "2018-12-01")

# convert to ts format
tax_ts <- as.ts(tax_tbl)

‍

Auch dies ist alles vom letzten Mal, ich habe es nur hier, damit Sie schnell loslegen können (einfach kopieren und einfügen)

# pretend we're in December 2017 and have to forecast the next twelve months
tax_ts_org <- window(tax_ts, end = c(2017, 12))

# estimate the required order of differencing
n_diffs <- nsdiffs(tax_ts_org)

# log transform and difference the data
tax_ts_trf <- tax_ts_org %>% 
  log() %>% 
  diff(n_diffs)

# embed the matrix
lag_order <- 6 # the desired number of lags (six months)
horizon <- 12 # the forecast horizon (twelve months)

tax_ts_mbd <- embed(tax_ts_trf, lag_order + 1) # embedding magic!

# do a train test split
y_train <- tax_ts_mbd[, 1] # the target
X_train <- tax_ts_mbd[, -1] # everything but the target

y_test <- window(tax_ts, start = c(2018, 1), end = c(2018, 12)) # the year 2018
X_test <- tax_ts_mbd[nrow(tax_ts_mbd), c(1:lag_order)] # the test set consisting
# of the six most recent values (we have six lags) of the training set. It's the
# same for all models.

‍

Cross-Validation und ihre (Fehl-)Anwendungen

Wie stimmt man die Hyperparameter eines ML-Modells ab, wenn man es mit Zeitreihendaten zu tun hat? Zuschauerhinweis: Das folgende Material könnte für manche verstörend sein.

Bevor wir die obige Frage beantworten, klären wir zunächst: Wie stimmt man die Hyperparameter eines ML-Modells ab, wenn die Daten keine Zeitdimension haben?

Im Allgemeinen gehen wir folgendermaßen vor:

1. Auswahl einiger Hyperparameter,
2. Festlegung eines sinnvollen Wertebereichs für jeden, und schließlich
3. Ermittlung der „besten“ Konfiguration, z. B. im Sinne der Minimierung eines Fehlermaßes.

Schritte 1) und 2) sind Trial-and-Error. In den meisten Fällen gibt es keinen analytischen Weg, um zu bestimmen, welche Hyperparameter in welcher Kombination am besten funktionieren. Zwar haben sich über die Jahre einige praktische Erfahrungswerte herausgebildet, im Grunde heißt es aber: ausprobieren und schauen, was bei deinen Daten funktioniert.

Und Schritt 3)? Die „beste“ Konfiguration findet man mithilfe verschiedener Resampling-Verfahren – eines der beliebtesten ist k-fold-Cross-Validation. Beliebt, weil k-fold-Cross-Validation recht verlässliche Schätzungen der Modellgüte liefert und die Daten effizient nutzt. Wie funktioniert das?

Zur Auffrischung: k-fold-Cross-Validation teilt die Daten in $k$ etwa gleich große Folds. Jeder Fold ist eine Teilmenge der Daten mit $N/k$ Beobachtungen. Das Modell wird dann auf $k-1$ Folds trainiert und auf dem $k$-ten, weggelassenen Fold validiert, indem man einen Loss (z. B. RMSE) berechnet. Dieser Prozess wird $k$-mal wiederholt, sodass jeder Fold einmal als Validierungsmenge dient.

Aus dieser Beschreibung ist sofort ersichtlich, warum das bei Zeitreihendaten problematisch ist. Wenn wir k-fold-Cross-Validation verwenden, trainieren wir das Modell zwangsläufig auf zukünftigen Daten, um die Vergangenheit vorherzusagen! Ich weiß nicht, wie es dir geht, aber das fühlt sich irgendwie nicht richtig an. Ganz zu schweigen davon, dass durch das zufällige Mischen der Daten vor dem Splitten die Zeitstruktur völlig zerstört wird.

Hold out – nicht so schnell

Was also tun? In der Zeitreihenökonometrie ist ein klassischer Ansatz: den letzten Teil der Zeitreihe als Holdout-Set reservieren. Das heißt, das Modell wird auf dem ersten Teil der Reihe trainiert und auf dem letzten Teil validiert. Im Wesentlichen ist das nichts anderes als eine Validierungsstrategie.

Dieser Ansatz ist einfach und intuitiv. Außerdem bewahrt er die zeitliche Reihenfolge der Daten. Aber er hat auch einen großen Nachteil: Er liefert keine robusten Schätzungen der Modellgüte, da wir das Modell nur ein einziges Mal testen – und das auf einem „Test“-Set, das willkürlich gewählt ist.

Ist das ein Problem?

Es kann eines sein: Zeitreihen repräsentieren vor allem Phänomene, die sich stetig verändern und weiterentwickeln. Insofern kann das Testset (also die „Zukunft“) systematisch anders sein als das Trainingsset (die „Vergangenheit“). Diese Veränderung im zugrunde liegenden Datenentstehungsprozess über die Zeit nennt man „Concept Drift“ und sie stellt eine ernsthafte Herausforderung in der Zeitreihenvorhersage dar.

Aber zurück zum Thema! Was wollen wir hier eigentlich erreichen? Wir wollen verlässliche Schätzungen der Modellgüte unseres ML-Modells bei Zeitreihendaten. Nur dann können wir mit einiger Sicherheit entscheiden, welche Hyperparameter wir wählen sollten. Doch der klassische Ansatz mit einem einzelnen Train-/Test-Split genügt dafür nicht.

Viel besser wäre es, mehrere Train-/Test-Splits zu verwenden. Im Kontext von Zeitreihen bedeutet das: Wir verschieben ein festes oder wachsendes Fenster über die Reihe, trainieren auf einem Abschnitt und sagen den nächsten voraus, berechnen den Fehler und verschieben das Fenster dann weiter (oder vergrößern es) und wiederholen den Prozess. Rob Hyndman hat dazu einen großartigen Beitrag geschrieben.

Dieser Ansatz, genannt Time Series Cross-Validation, ist effektiv – aber auch rechenintensiv. Stell dir vor, du hast 10 Hyperparameter-Kombinationen und testest jede mit 20 Train-/Test-Splits. Dann berechnest du insgesamt 200 Modelle. Je nach Modell und Datenmenge kann das einiges an Zeit in Anspruch nehmen.

Da wir beschäftigte Menschen in einer beschäftigten Welt sind, bleiben wir vorerst bei der Holdout-Strategie. Klingt einfach – aber wie macht man das konkret? Persönlich verwende ich dafür gerne die Funktion createTimeSlices aus dem Paket caret. Das ergibt vor allem Sinn, wenn du bereits mit einem caret-Workflow arbeitest oder viele unterschiedliche Serien modellierst. Wenn nicht, kannst du deine (Trainings-)Zeitreihe einfach so aufteilen, dass der letzte Abschnitt als Validierungsmenge dient.

Um zu sehen, wie die Funktion createTimeSlices funktioniert, führe sie zunächst eigenständig aus.

# we hold out the last 12 observations from our train set because this is also how far we want to predict into the future later
caret::createTimeSlices(
  1:nrow(X_train),
  initialWindow = nrow(X_train) - horizon,
  horizon = horizon,
  fixedWindow = TRUE
)

‍

Du siehst, dass unser Trainingssatz aufgeteilt wird in ein Trainings- und ein Validierungsset, basierend auf dem von uns angegebenen Prognosehorizont (12 Monate). Okay, richten wir das nun so ein, dass wir es in einem train-Workflow verwenden können:

tr_control <- caret::trainControl(
  method = 'timeslice',
  initialWindow = nrow(X_train) - horizon,
  fixedWindow = TRUE
)

‍

Mit trainControl als Grundlage richten wir als Nächstes ein Tuning-Grid ein. Auch wenn man hier sehr ausgeklügelt vorgehen kann, läuft es bei Random Forests meist auf zwei entscheidende Hyperparameter hinaus: die Anzahl der Bäume (ntree) und die Anzahl der Prädiktoren (mtry), die bei jedem Split im Baum zufällig ausgewählt werden.

Gute Werte für ntree liegen im Bereich einiger Hundert bis Tausend. Mehr Bäume können die Genauigkeit leicht verbessern, aber meistens nicht wesentlich. Da Random Forests kein hohes Overfitting-Risiko haben, hängt die Anzahl der verwendeten Bäume letztlich eher davon ab, wie viel Rechenleistung (oder Zeit) man zur Verfügung hat. Da das Tuning bei Zeitreihen typischerweise rechenintensiv ist, wählen wir hier 500 Bäume.

mtry hingegen ist in der Regel der spannendere Parameter. Er gibt an, wie viele Prädiktoren an jedem Knoten im Baum als Split-Kandidaten in Betracht gezogen werden. James et al. (2013) empfehlen als Richtwert $p/3$ oder $\sqrt{p}$ (wobei $p$ die Anzahl der Prädiktoren ist). Ich werde auch $p$ mit ins Grid aufnehmen – was dem Bagging entspricht (Nebenbemerkung: Wenn du diese Konzepte weiter vertiefen möchtest, sieh dir die Beiträge meiner Kolleg:innen zu Cross-Validation und Bagging an).

# caret actually only allows us to put one hyperparameter here
tune_grid <- expand.grid(
  mtry = c(
    ncol(X_train), # p
    ncol(X_train) / 3, # p / 3
    ceiling(sqrt(ncol(X_train))) # square root of p
  )
)

# let's see which hyperparameter the holdout method recommends
holdout_result <- caret::train(
  data.frame(X_train),
  y_train,
  method = 'rf',
  trControl = tr_control,
  tuneGrid = tune_grid
)

‍

Diese Methode empfiehlt mtry = 6. Also, sollten wir jetzt aufhören und unser Modell fitten? Noch nicht ganz!

Schon wieder Cross-Validation?

Erinnerst du dich an all die negativen Dinge, die wir über k-fold-Cross-Validation (CV) bei Zeitreihendaten gesagt haben? Technisch gesehen stimmt das alles immer noch. Aber jetzt kommt der überraschende Teil: Bergmeir et al. (2018) fanden heraus, dass k-fold-CV tatsächlich auf Zeitreihenmodelle angewendet werden kann – allerdings nur, wenn es sich um rein autoregressive Modelle handelt. Mit anderen Worten: Wir können k-fold-CV auf Zeitreihendaten anwenden, aber nur, wenn die Prädiktoren in unserem Modell zeitverzögerte Versionen der Zielvariable sind.

Wie das? Die Autoren stellen fest, dass k-fold-CV nur dann gültig ist, wenn die Fehler des Modells unkorreliert sind. Diese Bedingung ist erfüllt, wenn das Modell, das wir trainieren, die Daten gut abbildet. Anders gesagt: k-fold-CV ist gültig, wenn unser Modell ein gutes Modell „enthält“, d. h. ein Modell, dessen Gewichtungen eine Teilmenge der Gewichtungen unseres Modells sind. In diesem Fall ist k-fold-Cross-Validation sogar die bessere Wahl als die Holdout-Strategie!

Ich weiß nicht, wie es dir geht, aber für mich war das eine Überraschung. Da das eine äußerst nützliche Erkenntnis ist, lohnt sich die Wiederholung: Wenn die Residuen unseres (rein autoregressiven) Zeitreihenmodells seriell unkorreliert sind, dann kann und sollte k-fold-Cross-Validation der Holdout-Strategie vorgezogen werden.

Neugierig geworden? Dann wenden wir es doch direkt an – und zwar nicht nur einmal, sondern mehrfach!

tr_control <- trainControl(
  method = 'repeatedcv',
  number = 10, 
  repeats = 3
)

kfold_result <- caret::train(
  data.frame(X_train),
  y_train,
  method = 'rf',
  trControl = tr_control,
  tuneGrid = tune_grid
)

‍

Die Anwendung von k-fold-CV auf diese Zeitreihe schlägt einen Wert von mtry = 2 vor.

Interessant! Dann trainieren wir unseren Random Forest jetzt zweimal – einmal mit mtry = 2 und einmal mit mtry = 6. Danach vergleichen wir, welches Modell die bessere Vorhersage liefert!

mtrying it out!

# set up our empty forecast tibble
forecasts_rf <- tibble(
  mtry_holdout = rep(NA, horizon),
  mtry_kfold = rep(NA, horizon)
)

# collect the two mtry values from the tuning step
mtrys <- c(
  holdout_resultbestTune[[1]]
)

# train the model in a double loop
for (i in seq_len(length(mtrys))) {
  # start fresh for each mtry run
  y_train <- tax_ts_mbd[, 1] 
    X_train <- tax_ts_mbd[, -1] 

  # train the models
  for (j in seq_len(horizon)) {
    # set seed
    set.seed(2019)

    # fit the model
    fit_rf <- randomForest(X_train, y_train, mtry = mtrys[i])

    # predict using the test set
    forecasts_rf[j, i] <- predict(fit_rf, X_test)

    # here is where we repeatedly reshape the training data to reflect the time                            # distance corresponding to the current forecast horizon.
    y_train <- y_train[-1] 

    X_train <- X_train[-nrow(X_train), ] 
  }
}

‍

Auch diesmal müssen wir unsere Vorhersagen rücktransformieren, um die Genauigkeit auf dem Testset zu berechnen. Dafür verwende ich purrr.

last_observation <- as.vector(tail(tax_ts_org, 1))

forecasts <- forecasts_rf %>% 
  purrr::map_df(function(x) exp(cumsum(x)) * last_observation)

accuracies <- forecasts %>% 
  purrr::map(function(x) accuracy(x, as.vector(y_test))) %>%
  do.call(rbind, .)

‍

Und siehe da! k-fold-CV hat sich tatsächlich als besser erwiesen als unser Holdout-Ansatz. Sowohl RMSE als auch MAPE wurden reduziert. Wir haben sogar unser Ergebnis vom letzten Mal bestätigt, bei dem wir ebenfalls einen MAPE von 2,6 hatten. Zumindest in diesem Fall war es also völlig in Ordnung, Random Forest ohne spezielles Tuning einzusetzen.

‍

Visualisieren wir das Ganze auch noch:

plot <- tax_tbl %>% 
  filter(Date >= "2018-01-01") %>% 
  mutate(
    Ground_Truth = Value / 10000,
    Forecast_Holdout = forecastsmtry_kfold / 10000,
  ) %>% 
  ggplot(aes(x = Date)) +
  geom_line(aes(y = Value / 10000, linetype = "Truth")) +
  geom_line(aes(y = Forecast_Holdout, color = "Holdout")) +
  geom_line(aes(y = Forecast_Kfold, color = "k-fold CV")) +
  theme_minimal()+
  labs(
    title = "Forecast of the German Wage and Income Tax for the Year 2018",
    x = "Months",
    y = "Euros"
  ) +
  scale_color_manual(
    name = "Forecasts",
    values = c("Truth" = "black", "Holdout" = "darkblue", "k-fold CV" = "orange")
  ) +
  scale_linetype_manual(name = "Original", values = c("Truth" = "dashed")) +
  scale_x_date(date_labels = "%b %Y", date_breaks = "2 months")

‍

Wir sehen, dass die orangefarbene Linie, die die Vorhersagen des k-fold-CV-Modells repräsentiert, an mehreren Stellen enger an den echten Werten liegt.

Fazit (TL;DR)

Das Tuning von ML-Modellen auf Zeitreihendaten kann aufwendig sein – muss es aber nicht. Wenn das Modell ausschließlich endogene Prädiktoren verwendet, also Lags der Zielvariable, hast du Glück! Dann kannst du ganz beruhigt die bekannte und bewährte k-fold-Cross-Validation verwenden, um deine Hyperparameter abzustimmen. Wenn du tiefer einsteigen willst, lies dir unbedingt das zitierte Originalpapier durch. Oder du schnappst dir einfach eine Zeitreihe deiner Wahl und schaust, ob du dein Modell durch ein bisschen Tuning verbessern kannst. So oder so – es ist auf jeden Fall besser, als die Steuererklärung zu machen!

Quellen

• James, Gareth, et al. An Introduction to Statistical Learning. Vol. 112. New York: Springer, 2013.
• Bergmeir, Christoph, und José M. Benítez. “On the use of cross-validation for time series predictor evaluation.” Information Sciences 191 (2012): 192–213.