flights-delay-simplified

1 Hintergrund und Forschungsfrage

Wir untersuchen die Forschungsfrage Was sind Prädiktoren von Flugverspätungen. Wir nutzen dep_delay als AV (Zielvariable), also als die Variable, die wir vorhersagen wollen.

Dazu verwenden wir lineare Modelle als Modellierungsmethoden.

2 Ressourcen und Vertiefung

Dieser Post knüpft an den Post zur explorativen Datenanalyse der Flugverspätungen an (es gibt auch hier, Teil 1 und hier, Teil 2 ein Video zu diesem EDA-Post). Dieser Post ist eine längere (nicht-verkürzte) Version der vorliegenden Fallstudie.

3 Setup

3.1 Pakete laden

Wirklich wichtig sind nur tidymodels und tidyverse. Die restlichen Pakete werden nur am Rande benötigt. Man sollte auch nur die Pakete laden, die man für die Analyse benötigt.

library("tidymodels")  # Train- und Test-Sample aufteilen
library("tidyverse")  # data wrangling
library("conflicted")  # Name clashes finden
library("easystats")  # stats made easy
library("DataExplorer")  # Data Vis

3.2 Daten laden: Flights 2023

Aus Gründen der Datenökonomie nutzen wir eine kleinere Version des Datensatz flights. Wir nutzen nicht mehr die Daten aus dem 2013, sondern die neueren Daten aus dem Jahr 2023.

library(nycflights23)
data(flights)

set.seed(42)  # Reproduzierbarkeit
flights <- 
  flights |> 
  sample_n(size = 3e4)  # "3e4" heißt "3 Mal 10 hoch 4"

Achtung: flights ist recht groß; die Regressionsmodelle können leicht ein paar Hundert Megabyte groß werden. Das bringt u.U. auch einen modernen Computer irgendwann ins Straucheln. Daher verringern wir aus Gründen der Einfachheit den Datensatz mit einem Zufallssample. Man beachte, dass die Präzision der Ergebnisse höher ist, wenn man nicht mit einem Zufallssample, sondern dem gesamten Datensatz arbeitet.

4 flights2: Nicht benötigte Variablen entfernen und ID hinzufügen

flights2 <-
  flights %>% 
  select(-c(year)) %>%   # "year" ist konstant und daher nutzlos
  drop_na(dep_delay) %>% 
  mutate(id = row_number()) %>%  # ID als fortlaufende Nummer
  select(id, everything())  # id nach vorne ziehen

Die ID hilft uns, einzelne Flüge (Beobachtungen) zu identifizieren, z.B. wenn wir einen bestimmten Flug näher analysieren oder entfernen wolllen.

Achtung: Die ID-Variable sollte man nicht als Prädiktor verwenden, da sie keine Information birgt.

5 Heimliche AVs entfernen

dep_delay ist unsere AV (Zielvariable), die Variable also, die wir vorhersagen.

arr_delay ist eine Variable, die uns nicht zur Verfügung steht, wenn es um die Vorhersage neuer, noch nicht gestarteter Flüge geht. Es wäre daher falsch, unser Modell auf einer Variable zu gründen, auf die wir später (beim Vorhersagen) nicht zurückgreifen können. Also raus damit.

Ähnliches gilt für ein paar andere Variablen:

heimliche_avs <- c("dep_time", "arr_time", "arr_delay", "air_time")

flights2_ohne_heimliche_avs <-
  flights2 |> 
  select(-any_of(heimliche_avs))

any_of ist nur eine Tipphilfe. Man hätte auch direkt die Namen der zu entfernenden Spalten eingeben können, etwa select(-dep_time, ...).

Check:

names(flights2_ohne_heimliche_avs)

 [1] "id"             "month"          "day"            "sched_dep_time"
 [5] "dep_delay"      "sched_arr_time" "carrier"        "flight"        
 [9] "tailnum"        "origin"         "dest"           "distance"      
[13] "hour"           "minute"         "time_hour"     

Eine vielleicht gute Nachricht ist, dass Sie sich in der Prüfung nicht um diese Frage kümmern müssen.

6 Aufteilung in Train- und Testsample

Der Hintergrund zur Idee der Aufteilung in Train- und Test-Stichprobe kann z.b. hier oder hier, Kapitel 15, nachgelesen werden.

set.seed(42)  # Reproduzierbarkeit
flights_split <- initial_split(flights2_ohne_heimliche_avs, 
                               strata = dep_delay)

7 flights_train2, flights_test2

set.seed(42)  # Reproduzierbarkeit
flights_train2 <- training(flights_split)
flights_test2 <- testing(flights_split)

Die “wirkliche Welt” (was immer das ist) besorgt die Aufteilung von Train- und Test-Sample für Sie automatisch. Sagen wir, Sie arbeiten für die Flughafen-Aufsicht von New York. Dann haben Sie einen Erfahrungsschatz an Flügen aus der Vergangenheit in Ihrer Datenbank (Train-Sample). Einige Tages kommt Ihr Chef zu Ihnen und sagt: “Rechnen Sie mir mal die zu erwartende Verspätung der Flüge im nächsten Monat aus!”. Da heute nicht klar ist, wie die Verspätung der Flüge in der Zukunft (nächsten Monat) sein wird, stellen die Flüge des nächsten Monats das Test-Sample dar.

Übrigens: In der Prüfung besorgt das Aufteilen von Train- und Test-Sample netterweise Ihr Dozent…

Check:

names(flights_train2)

 [1] "id"             "month"          "day"            "sched_dep_time"
 [5] "dep_delay"      "sched_arr_time" "carrier"        "flight"        
 [9] "tailnum"        "origin"         "dest"           "distance"      
[13] "hour"           "minute"         "time_hour"     

8 lm0: Nullmodell

Eigentlich nicht nötig, das Nullmodell, primär aus didaktischen Gründen berechnet, um zu zeigen, dass in diesem Fall \(R^2\) wirklich gleich Null ist.

lm0 <- lm(dep_delay ~ 1, data = flights_train2)
model_parameters(lm0)  # model_parameters zeit die (geschätzten) Regressionsgewichte (Betas)

Parameter   | Coefficient |   SE |         95% CI | t(21918) |      p
---------------------------------------------------------------------
(Intercept) |       13.95 | 0.38 | [13.21, 14.69] |    36.91 | < .001

Wir könnten anstatt model_parameters auch parameters nutzen; das ist der gleiche Befehl.

Allerdings gibt es den Befehl parameters in zwei Paketen, es käme also zu einem “Name Clash”. Das umgehen wir, indem wir model_parameter nutzen, und nicht parameters.

9 lm1: origin

lm1 <- lm(dep_delay ~ origin, data = flights_train2)
model_parameters(lm1)  

Parameter    | Coefficient |   SE |         95% CI | t(21916) |      p
----------------------------------------------------------------------
(Intercept)  |       15.61 | 0.67 | [14.30, 16.92] |    23.40 | < .001
origin [JFK] |        0.19 | 0.95 | [-1.68,  2.05] |     0.20 | 0.845 
origin [LGA] |       -4.64 | 0.91 | [-6.42, -2.85] |    -5.09 | < .001

Man vergleiche:

flights_train2 %>% 
  drop_na(dep_delay) %>% 
  group_by(origin) %>% 
  summarise(delay_avg = mean(dep_delay)) %>% 
  mutate(delay_delta = delay_avg - delay_avg[1])

# A tibble: 3 × 3
  origin delay_avg delay_delta
  <chr>      <dbl>       <dbl>
1 EWR         15.6       0    
2 JFK         15.8       0.186
3 LGA         11.0      -4.64 

Der Mittelwertsvergleich und das Modell lm1 sind faktisch informationsgleich.

Aber leider ist es um die Modellgüte nicht so gut bestellt (eigentlich eher “grottenschlecht”):

r2(lm1)

# R2 for Linear Regression
       R2: 0.002
  adj. R2: 0.002

lm1 ist so schlecht, wir löschen es gleich wieder…

rm(lm1)

10 lm2: All in

# NICHT AUSFÜHREN
#lm2_all_in <- lm(dep_delay ~ ., data = flights_train2)

Modell lm2_all_in ist hier keine gute Idee, da nominale Prädiktoren in Indikatorvariablen umgewandelt werden. Hat ein nominaler Prädiktor sehr viele Stufen (wie hier), so resultieren sehr viele Indikatorvariablen, was dem Regressionsmodell Probleme bereiten kann (bei mir hängt sich R auf). Besser ist es in dem Fall, die Anzahl der Stufen von nominalskalierten Variablen vorab zu begrenzen.

Bei kleineren Datensätzen (weniger Variablen, weniger Fälle) lohnt es sich aber oft, das “All-in-Modell” auszuprobieren, als Referenzmaßstab für andere Modelle.

11 lm3: Verlorenes Modell

Dieses Modell ging verloren. Wo ist es hin? Wie geht es ihm? Leider gibt es keine Antwort in dieser Fallstudie; aber vielleicht in einer anderen…

Einigen Datensätzen ging es leider ebenso.¹

12 lm4: Alle metrischen Variablen

Was sind noch mal unsere metrischen Variablen:

flights_train4 <- 
flights_train2 %>% 
  select(where(is.numeric)) |> 
  select(-month, -day, -flight)   # diese Variablen sind nicht wirklich metrisch


flights_train4 %>% 
  names()

[1] "id"             "sched_dep_time" "dep_delay"      "sched_arr_time"
[5] "distance"       "hour"           "minute"        

Ok, jetzt eine Regression mit diesen Variablen (ober ohne die ID-Variable):

lm4 <- lm(dep_delay ~ ., 
          data = (flights_train4 |> select(-id)))

r2(lm4)

# R2 for Linear Regression
       R2: 0.016
  adj. R2: 0.016

Tja, das \(R^2\) hat einen nicht gerade um …

model_parameters(lm4)

Parameter      | Coefficient |       SE |          95% CI | t(21914) |      p
-----------------------------------------------------------------------------
(Intercept)    |       -7.41 |     1.42 | [-10.19, -4.63] |    -5.22 | < .001
sched dep time |        0.02 |     0.02 | [ -0.01,  0.06] |     1.26 | 0.208 
sched arr time |   -3.32e-03 | 1.10e-03 | [ -0.01,  0.00] |    -3.01 | 0.003 
distance       |    4.03e-03 | 5.36e-04 | [  0.00,  0.01] |     7.53 | < .001
hour           |       -0.79 |     1.94 | [ -4.59,  3.01] |    -0.41 | 0.683 

Die Distanz zum Ziel ist offenbar der interessanteste Prädiktor.

13 flights_train5: Fehlenden Werte ersetzen

flights_train4 |> 
  describe_distribution() |> 
  select(Variable, n_Missing)

Variable       | n_Missing
--------------------------
id             |         0
sched_dep_time |         0
dep_delay      |         0
sched_arr_time |         0
distance       |         0
hour           |         0
minute         |         0

Glücklicherweise haben wir keine fehlende Werte.

Nur rein zu Übungszwecken: Falls es fehlende Werte gibt, man könnte sie z.B. so mit dem Median ersetzen:

flights_train5 <-
  flights_train4 |> 
  mutate(distance = replace_na(distance, median(distance, na.rm = TRUE)))

14 lm6: Wie lm5, aber ohne fehlende Werte

lm6 <-lm(dep_delay ~ ., data = flights_train5 |> select(-id) )

r2(lm6)

# R2 for Linear Regression
       R2: 0.016
  adj. R2: 0.016

model_parameters(lm6)

Parameter      | Coefficient |       SE |          95% CI | t(21914) |      p
-----------------------------------------------------------------------------
(Intercept)    |       -7.41 |     1.42 | [-10.19, -4.63] |    -5.22 | < .001
sched dep time |        0.02 |     0.02 | [ -0.01,  0.06] |     1.26 | 0.208 
sched arr time |   -3.32e-03 | 1.10e-03 | [ -0.01,  0.00] |    -3.01 | 0.003 
distance       |    4.03e-03 | 5.36e-04 | [  0.00,  0.01] |     7.53 | < .001
hour           |       -0.79 |     1.94 | [ -4.59,  3.01] |    -0.41 | 0.683 

15 flights_train6: Extremwerte entfernen

flights_train5 |> 
  select(where(is.numeric), -id) |> 
  plot_density()

Es sieht so aus, als wäre distance deutlich rechtsschief mit einigen Ausreißern.

Man könnte auch Boxplots betrachten, die auch gut Extremwerte visualisieren.

Eine gängige Methode, mit Extremwerten umzugehen, ist, alle Datenpunkte, die im Boxplot als alleinstehende Punkte gezeigt werden, durch den Median zu ersetzen. Achtung: Diese Methode ist nicht perfekt! Es gibt viel sophistiziertere Methoden.

Wir ersetzen dabei alle Werte von distance, für die gilt, dass sie größer sind als Q3 + 1.5*IQR.

Q3:

flights_train5 |> 
  summarise(distance_q3 = quantile(distance, prob = .75))

# A tibble: 1 × 1
  distance_q3
        <dbl>
1        1183

IQR:

flights_train5 |> 
  summarise(distance_iqr = IQR(distance))

# A tibble: 1 × 1
  distance_iqr
         <dbl>
1          704

Der Grenzwert ist also:

(# Q3
  flights_train5 |> 
  summarise(distance_iqr = quantile(distance, prob = .75)) 
  ) +
  1.5 * 
 (# IQR
  flights_train5 |> 
  summarise(distance_iqr = IQR(distance)) 
 )

  distance_iqr
1         2239

Der Median von distance beträgt übrigens:

 flights_train5 |> 
  summarise(distance_md = median(distance)) 

# A tibble: 1 × 1
  distance_md
        <dbl>
1         762

flights_train6 <-
  flights_train5 |> 
  mutate(distance = 
           case_when(distance > 2239 ~ 762,
                     TRUE ~ distance))

Grob auf Deutsch übersetzt:

Wenn ein Flug eine distance von mehr als 326 Minuten hat, dann sei die airtime gleich 122, ansonsten immer (“TRUE”) ist airtime gleich distance, bleibt also, wie sie war.

16 lm7: Wie lm5, aber ohne Extremwerte für distance

lm7 <-lm(dep_delay ~ ., data = flights_train6)

r2(lm7)

# R2 for Linear Regression
       R2: 0.018
  adj. R2: 0.017

Tja….

17 R2 im Testsample

\(R^2\) kann man übrigens auch so berechnen:

Zuerst fügen wir die Vorhersagen zum Datensatz hinzu:

flights_train7_pred <- 
  flights_train6 %>% 
  mutate(lm7_pred = predict(lm7, newdata = flights_train6))  

Dann berechnen wir das R-Quadrat mit der Funktion rsq (wie “r squared”, R-Quadrat) anhand der beiden relevanten Spalten:

flights_train7_pred %>% 
  rsq(truth = dep_delay,
      estimate = lm7_pred)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rsq     standard      0.0176

Das ist nützlich, wenn man R-Quadrat auf Basis eigener Vorhersagen, im Test-Sample also, berechnen will.

Berechnen wir jetzt die Modellgüte im Testsample.

Fügen wir die Vorhersagewerte dem Testsample dazu:

flights_test4_pred <-
  flights_test2 %>% 
  mutate(pred_lm7 = predict(lm7, newdata = flights_test2))

Check:

flights_test4_pred %>% 
  select(id, dep_delay, pred_lm7) %>%
  head()

# A tibble: 6 × 3
     id dep_delay pred_lm7
  <int>     <dbl>    <dbl>
1    13         3     4.46
2    15        19    11.6 
3    28        -8    33.0 
4    30        -4    16.2 
5    31       160    20.8 
6    37        -8     2.79

flights_test4_pred |> 
  rsq(truth = dep_delay,
      estimate = pred_lm7)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rsq     standard      0.0184

Am schwierigsten ist es, bei den ganzen Nummerierungen nicht durcheinander zu kommen. Hier könnte es sich lohnen, ein übersichtlicheres Verfahren einzuführen (mit den Kosten höherer Komplexität).

18 Einreichen

Das beste Modell im Train-Sample reichen wir ein; in diesem Fall lm7.

submission_df <- 
flights_test4_pred |> 
  select(id, pred = pred_lm7)  # gleich umbenennen in "pred"

Check:

head(submission_df)

# A tibble: 6 × 2
     id  pred
  <int> <dbl>
1    13  4.46
2    15 11.6 
3    28 33.0 
4    30 16.2 
5    31 20.8 
6    37  2.79

Das sieht erstmal gut aus.

write_csv(submission_df, file = "Sauer_Sebastian_0123456_Prognose.csv")

19 Was noch?

Ein nächster Schritt könnte sein, sich folgende Punkte anzuschauen:

• Nominale Variablen berücksichtigen
• Variablen mit der höchsten Korrelation berücksichtigen
• Interaktionseffekte berücksichtigen

Eine Faustregel zu Interaktionen lautet: Wenn zwei Variablen jeweils einen starken Haupteffekt haben, lohnt es sich u.U., den Interaktionseffekt anzuschauen (vgl. Gelman & Hill, 2007, S. 69).

Footnotes

1. Hintergrund ist, dass diese Fallstudie eine vereinfachte und verkürzte Version einer ähnlichen Fallstudie ist.