1 Fallstudie NYC Flights

Aufgabe: Vorhersage von Verspätungen der Flüge von den Flughäfen von New York City im Jahr 2013.

1.1 Vorbereitung

Pakete laden:

library(mosaic)
library(tidyverse)
library(lubridate)
library(corrr)
library(caret)
library(doMC)
library(ranger)
library(sjmisc)

Daten laden:

library(nycflights13)
data(flights)
glimpse(flights)
#> Observations: 336,776
#> Variables: 19
#> $ year           <int> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2…
#> $ month          <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ day            <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ dep_time       <int> 517, 533, 542, 544, 554, 554, 555, 557, 557, 558,…
#> $ sched_dep_time <int> 515, 529, 540, 545, 600, 558, 600, 600, 600, 600,…
#> $ dep_delay      <dbl> 2, 4, 2, -1, -6, -4, -5, -3, -3, -2, -2, -2, -2, …
#> $ arr_time       <int> 830, 850, 923, 1004, 812, 740, 913, 709, 838, 753…
#> $ sched_arr_time <int> 819, 830, 850, 1022, 837, 728, 854, 723, 846, 745…
#> $ arr_delay      <dbl> 11, 20, 33, -18, -25, 12, 19, -14, -8, 8, -2, -3,…
#> $ carrier        <chr> "UA", "UA", "AA", "B6", "DL", "UA", "B6", "EV", "…
#> $ flight         <int> 1545, 1714, 1141, 725, 461, 1696, 507, 5708, 79, …
#> $ tailnum        <chr> "N14228", "N24211", "N619AA", "N804JB", "N668DN",…
#> $ origin         <chr> "EWR", "LGA", "JFK", "JFK", "LGA", "EWR", "EWR", …
#> $ dest           <chr> "IAH", "IAH", "MIA", "BQN", "ATL", "ORD", "FLL", …
#> $ air_time       <dbl> 227, 227, 160, 183, 116, 150, 158, 53, 140, 138, …
#> $ distance       <dbl> 1400, 1416, 1089, 1576, 762, 719, 1065, 229, 944,…
#> $ hour           <dbl> 5, 5, 5, 5, 6, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 6…
#> $ minute         <dbl> 15, 29, 40, 45, 0, 58, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
#> $ time_hour      <dttm> 2013-01-01 05:00:00, 2013-01-01 05:00:00, 2013-0…

1.2 Explorative Datenanalyse

1.2.1 Wie ist die Verspätung verteilt?

Es gibt zwei Variablen, die Verspätung anzeigen: arr_delay (Ankunft) und dep_delay.

favstats(arr_delay ~ 1, data = flights)
favstats(dep_delay ~ 1, data = flights)

Nehmen wir arr_delay, da die Streuung in dieser Variable höher ist.

1.2.2 VERTIEFUNG

Möchte man einen Befehl auf mehrere Spalten anwenden, so kann man dafür den Befehl map() verwendet. map() führt ein Befehl auf jede Spalte eines Dataframes aus. Damit man da Ergebnis in Form eines Dataframes (Tabelle) bekommt, fügt man _df an map() an:

flights %>% 
  select(arr_delay, dep_delay) %>% 
  map_df(favstats)

Für den IQR:

flights %>% 
  select(arr_delay, dep_delay) %>% 
  drop_na() %>% 
  map_df(iqr)

1.2.3 Visualisierung der Verteilung

gf_histogram( ~ arr_delay, data = flights)
gf_violin(arr_delay ~ 1, data = flights)

Aufgrund des langen rechten Randbereichs (hohe Verspätungswerte) ist das Diagramm nicht sher hilfreich.

Begrenzen wir uns besser auf den “inneren” Teil der Flüge (was die Verspätung betrifft).

flights %>% 
  filter(arr_delay < 120) %>% 
  gf_violin(arr_delay ~ 1, data = .)

1.2.4 Saisonale Effekte

Es gibt sehr viele potenzielle Ursachen für die Verspätung eines Flugzeugs bzw. eines Flugs. Zu einigen Kandidaten liegen uns Daten vor. Eine naheliegende (obwohl nicht tiefer theoretisch fundierte) Annahme ist, dass es saisonale Einflüsse auf die Verspätung gibt. So könnte Schnee im Winter oder Weihnachtsstress zum Jahreswechsel für Verspätung sorgen. Am Wochenende sind die Menschen entspannter und es wird weniger gereist. Daher könnte es Samstags und Sonntags zu weniger Verspätung kommen.

1.2.4.1 Nach Jahreszeiten

Berechnen wir die Jahreszeiten:

flights2 <- flights %>% 
  mutate(season = case_when(
    month %in% c(11, 12, 1, 2, 3) ~ "winter",
    month %in% c(6,7,9) ~ "summer",
    month %in% c(4, 5) ~ "spring",
    TRUE ~ "autumn"
  ))

Verspätungen nach Jahreszeiten:

favstats(arr_delay ~ season, data = flights2)

Im Sommer ist die Verspätung am höchsten. Vielleicht ist es besser, gleich auf Monate hin zu untersuchen:

favstats(arr_delay ~ month, data = flights2)

Tatsächlich ist die Verspätung im Mittelwert am höchsten im Juni und Juli. Dabei ist zu beachten, dass die mediane Verspätung nur im Dezember positiv ist: Nur im Dezember haben die Flüge in New York im Median eine Verspätung.

1.2.4.2 Weihnachten

Liegt es an Weihnachten? Schauen wir uns die Tage im Dezember (und Januar?) genauer an. Dazu berechnen wir zuerst einen Spalte, die den Tag (und die Woche) berechnet.

flights3 <- flights2 %>% 
  mutate(dayinyear = yday(time_hour),
         day_id = 365-(365-dayinyear),
         week = week(time_hour))
flights3 %>% 
  filter( (time_hour > "2013-11-30 23:59:59") ) %>% 
  group_by(dayinyear) %>% 
  summarise(arr_delay = median(arr_delay, na.rm = TRUE)) %>% 
  gf_line(arr_delay ~ dayinyear, data = .)

Etwa zwei Wochen vor Jahresende, also noch deutlich vor den Feiertagen, kommt es zu den Verspätungsspitzen. Ob zu dieser Zeit die meisten Menschen in den Weihnachtsurlaub fliegen? Insgesamt lässt diese Betrachtung offenbar keine starken Schlüsse zu.

1.2.4.3 Wochenende vs. Werktage

Vielleicht sind die Wochentage die entspannten Tage ohne Verspätung? Schauen wir nach. Man beachte, dass die Woche in Amerika mit Sonntag (1) beginnt, demzufolge ist der Samstag der 7. Tag.

flights3 %>% 
  mutate(weekend = if_else(wday(time_hour) %in% c(1,7), TRUE, FALSE)) %>% 
  group_by(weekend) %>% 
  summarise(arr_delay_md = median(arr_delay, na.rm = TRUE))

Aha, das sind 4 Minuten weniger im Median am Wochenende (im Vergleich zu werktags). Im Verhältnis zur Streuung von 31 Minuten (IQR) ist das nicht Nichts, aber auch nicht die Welt.

1.2.4.4 Verspätung pro Tag

flights4 <- flights3 %>% 
  mutate(day_rounded = round_date(time_hour, "day"))

flights4 %>% 
  group_by(day_rounded) %>% 
  summarise(arr_delay_md = median(arr_delay, na.rm = TRUE)) %>% 
  gf_line(arr_delay_md ~ day_rounded, data = .) %>% 
  gf_smooth()

Die Spitzen sind so nicht direkt erschließbar. Betrachten wir abschließend die Verspätungen pro Woche.

flights4 %>% 
  mutate(week_rounded = round_date(time_hour, "week")) %>% 
  group_by(week_rounded) %>% 
  summarise(arr_delay_md = median(arr_delay, na.rm = TRUE)) %>% 
  gf_line(arr_delay_md ~ week_rounded, data = .) %>% 
  gf_smooth()

Der Zacken im Juli könnte mit dem Nationalfeierag in den USA zusammenhängen. Lassen wir diese Untersuchungen an dieser Stelle.

1.2.5 Wetter

Die Wetterdaten sind in einer anderen Tabelle (weather), auch im Paket nycflights13 gespeichert. Über Datum/Zeit können wir die Wetterdaten mit den Flugdaten zusammenführen. Dabei begnügen wir uns mit einer tagesgenauen Präzision, da die Wetterdaten nicht jede Stunde (Minute, Sekunde) abdecken.

data(weather)
glimpse(weather)
#> Observations: 26,115
#> Variables: 15
#> $ origin     <chr> "EWR", "EWR", "EWR", "EWR", "EWR", "EWR", "EWR", "EWR…
#> $ year       <dbl> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013,…
#> $ month      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ day        <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ hour       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 15, 16, 17…
#> $ temp       <dbl> 39.02, 39.02, 39.02, 39.92, 39.02, 37.94, 39.02, 39.9…
#> $ dewp       <dbl> 26.06, 26.96, 28.04, 28.04, 28.04, 28.04, 28.04, 28.0…
#> $ humid      <dbl> 59.37, 61.63, 64.43, 62.21, 64.43, 67.21, 64.43, 62.2…
#> $ wind_dir   <dbl> 270, 250, 240, 250, 260, 240, 240, 250, 260, 260, 260…
#> $ wind_speed <dbl> 10.35702, 8.05546, 11.50780, 12.65858, 12.65858, 11.5…
#> $ wind_gust  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
#> $ precip     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
#> $ pressure   <dbl> 1012.0, 1012.3, 1012.5, 1012.2, 1011.9, 1012.4, 1012.…
#> $ visib      <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 1…
#> $ time_hour  <dttm> 2013-01-01 01:00:00, 2013-01-01 02:00:00, 2013-01-01…

weather2 <- weather %>% 
  mutate(date_time = round_date(time_hour, "hour"),
         day_rounded = round_date(time_hour, "day")) %>% 
  group_by(day_rounded) %>% 
  summarise_at(vars(temp, humid, wind_speed, precip, visib), median, na.rm = TRUE)
flights5 <- flights4 %>% 
  inner_join(weather2, by = c("day_rounded" = "day_rounded"))

Wie ist die Korrelation der Wetterdaten mit der Verspätung?

flights5 %>% 
  select(temp, humid, wind_speed, precip, visib, arr_delay) %>% 
  correlate() %>% 
  focus(arr_delay)

Gut, etwas Zusammenhang mit Luftfeuchtigkeit (humid), aber ansonsten nicht viel zu sehen. Apropos sehen: Schlechte Sicht geht mit weniger Verspätung einher (?).

1.3 Modellierung

1.3.1 Datensatz bereinigen

1.3.1.1 Variablen ohne Varianz

Variablen ohne Varianz sind wertlos für die Vorhersage, also entfernen wir sie. caret bietet dazu nearZeroVar. Natürlich kann man sich auch die Daten mit bloßem Auge ansehen, dann fällt auf, dass year den konstanten Wert 2013 aufweist.

flights6 <- flights5 %>% 
  select(-year)

1.3.1.2 Fehlende Werte

Probieren wir die rabiate Methode:

flights7 <- flights6 %>% 
  drop_na()

Wie viel Prozent der Fälle haben wir verloren?

nrow(flights7)/nrow(flights6)
#> [1] 0.9719892

Etwa 3%, das verschmerzen wir.

1.3.1.3 Z-Skalieren

Für viele Algorithmen ist es nötig (z.B. neuronale Netze), die Prädiktoren vorab zu standardisieren hinsichtlich Mittelwert und Streuung (z-Transformation). Das kann man z.B. so erreichen (via Paket sjmisc):

flights7a <- std(flights7, suffix = "")

1.3.2 Kreuzvalidierungsmethode

Um Überanpassung zu vermeiden, verwenden wir eine 5-fach-Kreuzvalidierung.

my_crossval <- trainControl(method = "cv",
                            number = 5,
                            allowParallel = TRUE,
                            verboseIter = FALSE)

allowParallel erlaubt die Verwendung mehrerer Rechenkerne, sofern initialisiert:

doMC::registerDoMC(cores = 2)

Achtung: Verdoppeln wir die Anzahl der Kerne, verdoppeln wir damit auch die Menge des benötigten Speichers.

1.3.3 Datensatz reduzieren

Große Datensätze bringen einen Rechner leicht aus der Ruhe. Besonders kategoriale Variablen mit vielen Stufen sind schwierig, da sie (manuell oder je nach Funktion automatisch) in Dummy-Variablen umgewandelt werden müssen.

Begrenzen wir uns daher auf metrische Variablen.

flights8 <- flights7a %>% 
  select_if(is.numeric)

Außerdem dürfen wir nicht vergessen, die andere Verspätungsvariable zu entferne (dep_delay).

flights9 <- flights8 %>% 
  select(-dep_delay)

1.3.4 Redundante Variablen

Haben wir noch redundante Variablen?

findLinearCombos(flights9)
#> $linearCombos
#> $linearCombos[[1]]
#> [1] 12  4 11
#> 
#> $linearCombos[[2]]
#> [1] 14 13
#> 
#> 
#> $remove
#> [1] 12 14

Ja!

Entfernen wir sie:

flights9a <- flights9 %>% 
  select(-c(12,14))

1.3.5 Datensatz aufteilen

Teilen wir den Datensatz zu 80% in einen Übungsteil bzw. zu 20% in einen Testdatensatz auf.

n_uebung <- round(.8 * nrow(flights9a), digits = 0)

uebung_index <- sample(1:nrow(flights9a), size = n_uebung)

uebung_df <- filter(flights9a, row_number() %in% uebung_index)
test_df <- filter(flights9a, !(row_number() %in% uebung_index))

Die Gesamtfallzahl muss der Summe aus Übungs- und Test-Datensatz entsprechen:

(nrow(uebung_df) + nrow(test_df)) == nrow(flights9a)
#> [1] TRUE

Passt.

1.3.6 Modell 1 - Regression

Beginnen wir mit einer einfachen Regression (ohne Interaktionen, Polynome, etc.).1

Eine Regression hat keine Tuningparameter.

start <- Sys.time()
lm_fit1 <- train(arr_delay ~ .,
                 data = uebung_df,
                 method = "lm",
                 trControl = my_crossval)
end <- Sys.time()

(time_taken_lm1 <- end - start)
#> Time difference of 16.67416 secs


#saveRDS(lm_fit1, file = "lm_fit1.rds")

Ohne die Begrenzung auf numerische Variablen hat meine Maschine (16GB Speicher) einen Asthmaanfall bekommen und ist steckengeblieben.2

Das erzeugte Modell hatte in der Datei eine Größe von ca. 360MB.

Die Koeffizienten des Modells lassen sich auf übliche Weise bestimmen:

summary(lm_fit1)
#> 
#> Call:
#> lm(formula = .outcome ~ ., data = dat)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2.1075 -0.4652 -0.1807  0.1795 25.6674 
#> 
#> Coefficients:
#>                  Estimate Std. Error  t value Pr(>|t|)    
#> (Intercept)    -1.736e-06  1.761e-03   -0.001  0.99921    
#> month           8.164e+00  3.014e-01   27.082  < 2e-16 ***
#> day             7.014e-01  2.547e-02   27.542  < 2e-16 ***
#> dep_time        7.888e-01  6.099e-03  129.345  < 2e-16 ***
#> sched_dep_time -5.888e-01  4.338e-02  -13.574  < 2e-16 ***
#> arr_time       -2.634e-01  2.903e-03  -90.755  < 2e-16 ***
#> sched_arr_time  8.686e-02  3.582e-03   24.252  < 2e-16 ***
#> flight          3.899e-02  2.017e-03   19.328  < 2e-16 ***
#> air_time        1.643e+00  1.410e-02  116.532  < 2e-16 ***
#> distance       -1.664e+00  1.415e-02 -117.593  < 2e-16 ***
#> hour            1.158e-01  4.282e-02    2.705  0.00683 ** 
#> dayinyear      -8.561e+00  3.172e-01  -26.992  < 2e-16 ***
#> week            3.548e-01  9.062e-02    3.915 9.03e-05 ***
#> temp            7.449e-02  2.372e-03   31.405  < 2e-16 ***
#> humid           1.530e-01  2.201e-03   69.524  < 2e-16 ***
#> wind_speed      9.915e-02  1.958e-03   50.640  < 2e-16 ***
#> precip         -2.175e-03  1.933e-03   -1.125  0.26042    
#> visib           6.001e-03  2.198e-03    2.730  0.00632 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.9006 on 261512 degrees of freedom
#> Multiple R-squared:  0.1871, Adjusted R-squared:  0.1871 
#> F-statistic:  3541 on 17 and 261512 DF,  p-value: < 2.2e-16

Die Prädiktorenrelevanz kann man über varImp() abfragen.

varImp(lm_fit1)
#> lm variable importance
#> 
#>                Overall
#> dep_time       100.000
#> distance        90.834
#> air_time        90.007
#> arr_time        69.903
#> humid           53.345
#> wind_speed      38.617
#> temp            23.616
#> day             20.602
#> month           20.244
#> dayinyear       20.173
#> sched_arr_time  18.036
#> flight          14.196
#> sched_dep_time   9.709
#> week             2.176
#> visib            1.252
#> hour             1.232
#> precip           0.000

1.3.7 Modell 2 - Random Forest

Im Gegensatz zur Regression gibt es bei Random-Forest-Modellen Tuningparameter, und zwar die Anzahl der Variablen pro Baum, hier mit .mtry bezeichnet. Eine Faustregel für diesen Parameter ist \(\sqrt(k)\), hier also etwa oder 6.

rf_grid <- data.frame(
  .mtry = c(4, 5, 6, 7),
  .splitrule = "variance",
  .min.node.size = 5)

rf_grid

Um Zeit zu sparen, verringern wir die Stichprobengröße auf 1000:

uebung_df_small <- sample_n(uebung_df, size = 1000)

Dann berechnen wir das Modell; gibt man keine Hinweise auf Variation von Tuningparametern, so wählt die Funktion Standardwerte.

start <- Sys.time()
rf_fit1 <- train(arr_delay ~ .,
                 data = uebung_df_small,
                 method = "ranger",
                 trControl = my_crossval)
end <- Sys.time()

(time_taken <- end - start)
#> Time difference of 17.81419 secs

#saveRDS(rf_fit1, file = "lm_fit1.rds")
#readRDS("lm_fit1.rds")

Einen Überblick über das berechnete Modell kann man sich so ausgeben lassen:

rf_fit1
#> Random Forest 
#> 
#> 1000 samples
#>   17 predictor
#> 
#> No pre-processing
#> Resampling: Cross-Validated (5 fold) 
#> Summary of sample sizes: 800, 800, 799, 801, 800 
#> Resampling results across tuning parameters:
#> 
#>   mtry  splitrule   RMSE       Rsquared   MAE      
#>    2    variance    0.7633726  0.3607519  0.5009789
#>    2    extratrees  0.7907793  0.3392968  0.5121332
#>    9    variance    0.6610500  0.5292866  0.4392290
#>    9    extratrees  0.6827074  0.5194783  0.4581119
#>   17    variance    0.6429513  0.5607234  0.4243621
#>   17    extratrees  0.6495064  0.5659727  0.4361470
#> 
#> Tuning parameter 'min.node.size' was held constant at a value of 5
#> RMSE was used to select the optimal model using the smallest value.
#> The final values used for the model were mtry = 17, splitrule =
#>  variance and min.node.size = 5.

Im resultierenden Objekt sind eine Vielzahl von Informationen zu finden. So kann man sich den Modellkandidaten mit den besten Werten ausgeben lassen:

rf_fit1$bestTune

Aber was ist das Kriteriem, das optimiert wird?

rf_fit1$metric
#> [1] "RMSE"

Es wird nach dem Root Mean Square Error optimiert.

Weitere Infos zum Algorithmus bekommt man z.B. so:

modelLookup("ranger")

Wir sehen, dass das Modell drei Tuningparameter hat (wobei caret den Parameter min.node.size konstant hielt).

plot(rf_fit1)

Möchte man eine bestimmte Anzahl an Kandidatenmodelle prüfen lassen, so kann man das mit tuneLength tun:

start <- Sys.time()
rf_fit2 <- train(arr_delay ~ .,
                 data = uebung_df_small,
                 method = "ranger",
                 trControl = my_crossval,
                 tuneLength = 4)
end <- Sys.time()


(time_taken <- end - start)
#> Time difference of 22.06099 secs

saveRDS(rf_fit2, file = "lm_fit2.rds")
rf_fit2
#> Random Forest 
#> 
#> 1000 samples
#>   17 predictor
#> 
#> No pre-processing
#> Resampling: Cross-Validated (5 fold) 
#> Summary of sample sizes: 800, 801, 801, 799, 799 
#> Resampling results across tuning parameters:
#> 
#>   mtry  splitrule   RMSE       Rsquared   MAE      
#>    2    variance    0.7623641  0.3546775  0.5016084
#>    2    extratrees  0.7818499  0.3589807  0.5089059
#>    7    variance    0.6833994  0.4963581  0.4493724
#>    7    extratrees  0.6949830  0.5088313  0.4633939
#>   12    variance    0.6557907  0.5389356  0.4313334
#>   12    extratrees  0.6622967  0.5540862  0.4441157
#>   17    variance    0.6490689  0.5545829  0.4246735
#>   17    extratrees  0.6419074  0.5859532  0.4320974
#> 
#> Tuning parameter 'min.node.size' was held constant at a value of 5
#> RMSE was used to select the optimal model using the smallest value.
#> The final values used for the model were mtry = 17, splitrule =
#>  extratrees and min.node.size = 5.

Mit tuneGrid kann man die Werte der Modellkandidaten genau einstellen.

start <- Sys.time()
rf_fit3 <- train(arr_delay ~ .,
                 data = uebung_df_small,
                 method = "ranger",
                 trControl = my_crossval,
                 tuneGrid = rf_grid)
end <- Sys.time()


(time_taken <- end - start)
#> Time difference of 8.216766 secs

# saveRDS(rf_fit3, file = "lm_fit3.rds")
rf_fit3
#> Random Forest 
#> 
#> 1000 samples
#>   17 predictor
#> 
#> No pre-processing
#> Resampling: Cross-Validated (5 fold) 
#> Summary of sample sizes: 799, 800, 800, 800, 801 
#> Resampling results across tuning parameters:
#> 
#>   mtry  RMSE       Rsquared   MAE      
#>   4     0.7116300  0.4554067  0.4709582
#>   5     0.6922898  0.4903751  0.4593296
#>   6     0.6816250  0.5020712  0.4519450
#>   7     0.6710419  0.5200150  0.4457520
#> 
#> Tuning parameter 'splitrule' was held constant at a value of
#>  variance
#> Tuning parameter 'min.node.size' was held constant at a value
#>  of 5
#> RMSE was used to select the optimal model using the smallest value.
#> The final values used for the model were mtry = 7, splitrule =
#>  variance and min.node.size = 5.

1.3.8 Modell 3 - Neuronales Netz

Neuronale Netze benötigen große Datensätze, daher ist unser kleiner Datensatz mit \(n=1000\) sicher zu klein (gerade in Anbetracht zur Zahl der Features). Aus Vergleichsbarkeitsgründen und um die Rechenkosten zu schätzen, bietet es sich aber an, zunächst mit einem kleinen Datensatz zu arbeiten:

start <- Sys.time()
nn_fit1 <- train(arr_delay ~ .,
                 data = uebung_df_small,
                 method = "nnet",
                 trControl = my_crossval,
                 linout = TRUE)
#> # weights:  96
#> initial  value 1251.713470 
#> iter  10 value 752.392005
#> iter  20 value 629.992664
#> iter  30 value 510.727200
#> iter  40 value 370.835270
#> iter  50 value 303.223878
#> iter  60 value 262.263392
#> iter  70 value 231.974569
#> iter  80 value 191.799761
#> iter  90 value 163.597567
#> iter 100 value 143.767454
#> final  value 143.767454 
#> stopped after 100 iterations
end <- Sys.time()


(time_taken <- end - start)
#> Time difference of 3.906937 secs

saveRDS(nn_fit1, file = "nn_fit1.rds")

Der Parameter linout = TRUE verhindert eine Aktivierungsfunktion, die den Wertebereich auf [0,1] beschränken würde.

Das ging schnell. Vergrößen wir den Datensatz:

start <- Sys.time()
nn_fit2 <- train(arr_delay ~ .,
                 data = uebung_df,
                 method = "nnet",
                 trControl = my_crossval,
                 linout = TRUE)
#> # weights:  96
#> initial  value 296093.851805 
#> iter  10 value 206345.983463
#> iter  20 value 123079.214569
#> iter  30 value 100052.275388
#> iter  40 value 90280.334351
#> iter  50 value 81869.501514
#> iter  60 value 73615.105786
#> iter  70 value 62928.199709
#> iter  80 value 52795.472294
#> iter  90 value 47514.499738
#> iter 100 value 43863.841483
#> final  value 43863.841483 
#> stopped after 100 iterations
end <- Sys.time()


(time_taken <- end - start)
#> Time difference of 12.04127 mins

# saveRDS(nn_fit2, file = "nn_fit2.rds")

Es gibt verschiedene Implementierungen von neuronalen Netzen, die in caret angesteuert werden können, z.B. neuralnet. Es verfügt über 3 Modellparameter:

modelLookup("neuralnet")
getModelInfo("neuralnet")
#> $neuralnet
#> $neuralnet$label
#> [1] "Neural Network"
#> 
#> $neuralnet$library
#> [1] "neuralnet"
#> 
#> $neuralnet$loop
#> NULL
#> 
#> $neuralnet$type
#> [1] "Regression"
#> 
#> $neuralnet$parameters
#>   parameter   class                    label
#> 1    layer1 numeric #Hidden Units in Layer 1
#> 2    layer2 numeric #Hidden Units in Layer 2
#> 3    layer3 numeric #Hidden Units in Layer 3
#> 
#> $neuralnet$grid
#> function(x, y, len = NULL, search = "grid") {
#>                     if(search == "grid") {
#>                       out <- expand.grid(layer1 = ((1:len) * 2) - 1, layer2 = 0, layer3 = 0)
#>                     } else {
#>                       out <- data.frame(layer1 = sample(2:20, replace = TRUE, size = len),
#>                                         layer2 = sample(c(0, 2:20), replace = TRUE, size = len),
#>                                         layer3 = sample(c(0, 2:20), replace = TRUE, size = len))
#>                     }
#>                     out
#>                   }
#> 
#> $neuralnet$fit
#> function(x, y, wts, param, lev, last, classProbs, ...) {
#>                     colNames <- colnames(x)
#>                     dat <- if(is.data.frame(x)) x else as.data.frame(x)
#>                     dat$.outcome <- y
#>                     form <- as.formula(paste(".outcome ~",paste(colNames, collapse = "+")))
#>                     if(param$layer1 == 0) stop("the first layer must have at least one hidden unit")
#>                     if(param$layer2 == 0 & param$layer2 > 0) stop("the second layer must have at least one hidden unit if a third layer is specified")
#>                     nodes <- c(param$layer1)
#>                     if(param$layer2 > 0) {
#>                       nodes <- c(nodes, param$layer2)
#>                       if(param$layer3 > 0) nodes <- c(nodes, param$layer3)
#>                     }
#>                     neuralnet::neuralnet(form, data = dat, hidden = nodes, ...)
#>                   }
#> 
#> $neuralnet$predict
#> function(modelFit, newdata, submodels = NULL) {
#>                     newdata <- newdata[, modelFit$model.list$variables, drop = FALSE]
#>                     neuralnet::compute(modelFit, covariate = newdata)$net.result[,1]
#>                   }
#> 
#> $neuralnet$prob
#> NULL
#> 
#> $neuralnet$tags
#> [1] "Neural Network"
#> 
#> $neuralnet$sort
#> function(x) x[order(x$layer1, x$layer2, x$layer3),]

1.4 Vergleich der Modellgüten

Wie gut sagen die Modelle den Test-Datensatz vorher? Vergleichen wir die Modelle.

1.4.1 Prognosen für den Test-Datensatz berechnen

Erstellen wir uns einen Datensatz mit den Vorhersagen (und den beobachteten Werten):


test_preds <- test_df %>% 
  select(arr_delay) %>% 
  mutate(lm1_pred = predict(lm_fit1, newdata = test_df))

test_preds <- test_df %>% 
  select(arr_delay) %>% 
  mutate(lm1_pred = predict(lm_fit1, newdata = test_df),
         rf1_pred = predict(rf_fit2, newdata = test_df),
         rf2_pred = predict(rf_fit3, newdata = test_df),
         nn1_pred = predict(nn_fit1, newdata = test_df),
         nn2_pred = predict(nn_fit2, newdata = test_df))

Jetzt lassen wir uns typische Kennzahlen der Modellgüte ausgeben:

postResample(pred = test_preds$lm1_pred, obs = test_preds$arr_delay)
#>      RMSE  Rsquared       MAE 
#> 0.9037563 0.1906945 0.5425576

Der Vektor der Modellnamen lautet:

model_names <- names(test_preds)

Das wiederholen wir in einer Schleife für jedes Modell:

test_pred_df <- test_preds %>% 
  map_df(~ postResample(pred = ., obs = test_preds$arr_delay)) %>% 
  mutate(statistic = c("RMSE", "Rsquared", "MAE")) %>% 
  select(statistic, everything())

test_pred_df

Formen wir diese Tabelle in Langform (Normalform um):

test_pred_df_t <- test_pred_df %>% 
  gather(key = "model_name", value = "value", -c(statistic))

Eine andere Form wäre:

test_pred_df %>% 
  gather(key = "model_name", value = "value", -c(statistic)) %>% 
  spread(key = statistic, value = value)

1.4.2 Bestes Modell identifizieren

Der kleinste RMSE-Wert (nach dem Modell, dass als vorhergesagten Werte die beobachteten nimmt, also einen RSME von Null hat):

test_pred_df_t %>% 
  filter(statistic == "RMSE") %>% 
  top_n(2, wt = -value)

Der größte R^2-Wert:

test_pred_df_t %>% 
  filter(statistic == "Rsquared") %>% 
  top_n(2, wt = value)

1.4.3 Visualisieren

test_pred_df_t %>% 
  group_by(statistic) %>% 
  mutate(is_max = value == max(value),
         is_min = value == min(value)) %>% 
  ggplot(aes(y = model_name, x = value, color = is_max, shape = is_min)) +
  geom_point(size = 5) +
  facet_wrap(~statistic, scales = "free_x") +
  theme(legend.position = "bottom")

Damit hat das neuronale Netz “gewonnen”.

2 Fazit

Es darf nicht vergessen werden, dass wir nur einen Teil des Datensatzes verwendet haben - schlicht aus Gründen der komputationalen Kostensparung. Insofern sind die Modellgüten nur bedingt für bare Münze zu nehmen. Diese Fallstudie hat nur einen Teil der Möglichkeiten einer ernsthaften Modellierung aufgenommen, so dass die Ergebnisse schon aus diesem Grund mit einem großen Gramm Salz zu betrachten sind. Außerdem fanden nur relativ weniger Modell Eingang; es bleibt also offen, ob nicht andere Modelle “besser” sind. Beim Wort “besser” muss man immer im Kopf behalten, dass “besser” eine bedingte Aussage ist: besser vor dem Hintergrund gewisser anderer Modelle, gewisser Transformationen, gewisser Implementierungen, gewisser Stichprobenmerkmale, gewisser Implementierungsspezifika und so weiter.

  1. vgl. https://topepo.github.io/caret/train-models-by-tag.html#linear-regression↩︎

  2. vgl. https://stackoverflow.com/questions/51248293/error-vector-memory-exhausted-limit-reached-r-3-5-0-macos?rq=1↩︎