1 Pakete laden

library(tidyverse)  # data wrangling
library(fastDummies)  # nur für "Dummyisierung"
library(skimr)  # viele Statistiken auf einmal
library(corrr)  # komfortabel Korrelationen ausrechnen

2 Hintergrund und Ziel

Dieser Post zeigt einige mögliche/typische Schritte der explorativen Datenanalyse (EDA) im Hinblick auf die Forschungsfrage “Welche Variablen steht in Zusammenhang mit Flugverspätungen?”.

3 Daten laden

Der Datensatz kann z.B. hier bezogen werden:

flights <- read_csv("https://vincentarelbundock.github.io/Rdatasets/csv/nycflights13/flights.csv")

Der Datensatz ist relativ groß:

dim(flights)
#> [1] 336776     20

Alternativ findet man den Datensatz auch im Paket nycflights13.

data(flights, package = "nycflights13")

4 Was ist Verspätung?

Schauen wir uns den Datensatz mal näher an, um die Zielvariable “Verspätung” zu beleuchten.

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

Es gibt zwei Spalten, die auf Verspätung hinzielen:

flights %>% 
  select(contains("delay")) %>% 
  head()

Ein Blick in die Hilfe lässt uns mehr lernen, was die Spalten bedeuten: help(flights) (wenn das Paket nycflights13 geladen ist; alternativ kann man z.B. hier nachlesen):

dep_delay, arr_delay

Departure and arrival delays, in minutes. 
Negative times represent early departures/arrivals.

flights %>% 
  sample_n(size = 1000) %>% 
  ggplot() +
  aes(y = dep_delay, x = arr_delay) +
  geom_point()

flights %>%
  drop_na(dep_delay, arr_delay) %>% 
  summarise(sd(dep_delay),
            sd(arr_delay))

flights %>% 
  select(contains("delay")) %>% 
  drop_na() %>% 
  mutate(dep_delay = scale(dep_delay),  # z-Transformation
         arr_delay = scale(arr_delay)) %>%   # z-Transformation
  ggplot() +
  aes(x = arr_delay, y = dep_delay) +
  geom_bin2d() +
  geom_abline(linetype = "dashed",
              color = "grey60")

flights %>% 
  drop_na(dep_delay, arr_delay) %>% 
  summarise(cor(dep_delay, arr_delay))

5 Verteilung der Verspätung

flights %>% 
  ggplot() +
  aes(x = dep_delay) %>% 
  geom_density()

Ein sehr langer rechter Rand; die meisten Flüge sind nicht/kaum verspätet; aber einige wenige sind sehr stark verspätet.

Zentrale deskriptive Statisitken könnte man sich mit summary ausgeben lassen:

flights %>%
  filter(!is.na(dep_delay)) %>%  # keine fehlenden Werte
  summarise(depdelay_mean = mean(dep_delay),
            depdelay_sd = sd(dep_delay),
            depdelay_md = median(dep_delay),
            depdelay_iqr = IQR(dep_delay)) 

Oder man benutzt den Befehl across, der es erlaubt, eine oder mehrere Funktionen auf eine oder mehrere Spalten anzuwenden. In diesem Beispiel wenden wir mehrere Funktionen (adressiert mit .fns) auf eine Spalte (dep_delay), adressiert mit dem Argument .cols an. Außerdem kann man die Namen der resultierenden Spalten bestimmen mit dem Argument .names. In der geschweiften Klammer steht eine interne Variable, die den Namen der jeweils berechneten Funktion ({fn}) an den Namen der neu erstellten Spalte anfügt; in der Ausgabe sieht man das gut.

flights %>% 
  summarise(across(
    .cols = dep_delay,
    .fns = list(mean = mean, 
                md = median, 
                sd = sd, 
                iqr = IQR), na.rm = TRUE,
    .names = "depdelay_{fn}"
  ))

flights %>% 
  summarise(q75 = quantile(dep_delay, 
                           prob = .75, 
                           na.rm = TRUE))

flights %>% 
  summarise(depdelay_iqr = IQR(dep_delay, na.rm = TRUE))

grenzwert <- 11 + 1.5*16

flights2 <-
  flights %>% 
  mutate(is_extreme = case_when(
    dep_delay > 11 + 1.5 * 16 ~ TRUE, # Verspätung > 35 Min.
    dep_delay <= 35 ~ FALSE  # in den anderen Fällen (<= 35 Min.), dann kein Extremwert
  ))

flights %>% 
  summarise(sum(is.na(dep_delay)))  # fehlende Werte zählen

flights %>% 
  summarise(across(everything(), ~ sum(is.na(.x))))

flights %>% 
  drop_na(dep_delay) %>% 
  nrow()
#> [1] 328521

flights %>% 
  drop_na() %>% 
  nrow()
#> [1] 327346

flights3 <-
  flights2 %>% 
  drop_na() %>% 
  select(-year)

6 Deskriptive Statistiken

flights2 %>% 
  summarise(mean(dep_delay),
            sd(dep_delay),
            mean(arr_delay),
            sd(arr_delay))  # und so weiter

flights %>% 
  skim()

7 Korrelate von Verspätung

Schauen wir, welche Variablen mit dep_delay, der Verspätung der Flüge also, korrelieren.

flights3 %>% 
  select(where(is.numeric)) %>% 
  summarise(cor_month = cor(dep_delay, month),
            cor_day = cor(dep_delay, day),
            cor_dep_time = cor(dep_delay, dep_time))  # etc

flights3 %>% 
  select(where(is.numeric)) %>%  # nur die numerischen Spalten auswählen
  summarise(across(
    .cols = everything(),
    .fns = ~ cor(., dep_delay))) %>% 
  pivot_longer(everything()) %>%  # von breit auf lang
  arrange(-value)  # absteigend sortieren

flights2 %>% 
  select(where(is.numeric)) %>%  # alle metrischen Variablen, aber sonst keine
  correlate() %>%   # korreliere, was du hast
  focus(dep_delay) %>%  # beschränke (fokussiere) dich auf `dep_delay`
  arrange(-dep_delay)

flights2 %>% 
  select(where(negate(is.numeric))) %>% 
  names()
#> [1] "carrier"    "tailnum"    "origin"     "dest"       "time_hour" 
#> [6] "is_extreme"

flights2 %>% 
  select(where(negate(is.numeric))) %>% 
  summarise(across(
    .cols = everything(),
    .fns = list(ndist = n_distinct)
  )) %>% 
  pivot_longer(everything())

flights2 %>% 
  select(carrier) %>% 
  count(carrier, sort = TRUE) 

plot1 <- flights2 %>% 
  select(carrier) %>% 
  count(carrier, sort = TRUE) %>% 
  ggplot() +
  aes(y = carrier, x = n) +
  geom_point(color = "red") +
  geom_line(group = 1)

plot1

plot1 +
  scale_x_continuous(breaks = c(0, 25000, 50000, 75000),
                     limits = c(0, 100000),
                     labels = c("keine", "wenig", "mittel", "viel"))

flights2 <-
  flights2 %>% 
  mutate(carrier = factor(carrier)) %>%   # nicht character, sondern factor
  mutate(carrier_lump = fct_lump(carrier, n = 8))

flights2 %>% 
  count(carrier_lump)

flights2 %>% 
  select(carrier) %>% 
  count(carrier, sort = TRUE) %>% 
  ggplot() +
  aes(y = fct_reorder(carrier, n), x = n) +
  geom_point(color = "red") +
  geom_line(group = 1)

flights2 %>% 
  filter(!is_extreme) %>%  # keine Extremwerte
  ggplot() +
  aes(x = fct_reorder(carrier_lump, dep_delay, 
                      na.rm = TRUE), 
      y = dep_delay) %>% 
  geom_boxplot()

flights2 %>% 
  filter(!is_extreme) %>% 
  group_by(carrier_lump) %>% 
  summarise(dep_delay = median(dep_delay, na.rm = TRUE)) %>% 
  arrange(dep_delay)

flights2 <- 
  flights2 %>% 
  mutate(
    originJFK = case_when(
      origin == "JFK" ~ 1,
      origin != "JFK" ~ 0
    ),
    originLGA = case_when(
      origin == "LGA" ~ 1,
      TRUE ~ 0,
    )
  )

flights2 %>% 
  select(origin, originJFK, originLGA) %>% 
  slice(1:5)

flights2 %>% 
  select(origin, dep_delay) %>% 
  dummy_cols() %>%  # aus dem Paket `fastDummies`
  head()  # slice(1:6)

flights2 %>% 
  select(dep_delay, carrier) %>% 
  dummy_cols() %>%   # "Dummysierung"
  select(-carrier) %>% 
  pivot_longer(-dep_delay,
               names_to = "carrier",
               values_to = "value") %>% 
  group_by(carrier) %>% 
  summarise(cor_depdelay_carrier = cor(dep_delay, value,
                                       use = "complete.obs")) %>% 
  arrange(-abs(cor_depdelay_carrier)) %>% 
  filter(abs(cor_depdelay_carrier) > 0.10)

flights2 %>% 
  filter(!is_extreme) %>% 
  select(dep_delay, hour) %>% 
  mutate(hour = factor(hour)) %>% 
  ggplot() +
  aes(x = fct_reorder(hour, dep_delay,
                      na.rm = TRUE), 
      y = dep_delay) +
  geom_boxplot()

flights2 %>% 
  filter(!is_extreme) %>% 
  select(origin, dep_delay) %>% 
  ggplot() +
  aes(x = origin, y = dep_delay) %>% 
  geom_boxplot()

flights2 %>% 
  filter(!is_extreme) %>% 
  select(origin, dep_delay, hour) %>% 
  mutate(hour = factor(hour, levels = 5:23)) %>% 
  ggplot() +
  aes(x = hour, y = dep_delay) +
  geom_boxplot() +
  facet_wrap(~ origin)

flights3 %>% 
  select(carrier, tailnum, origin, dest, time_hour) %>% 
  slice(1:3)

flights3 %>% 
  select(where(~ !is.numeric(.))) %>% 
  names()
#> [1] "carrier"    "tailnum"    "origin"     "dest"       "time_hour" 
#> [6] "is_extreme"

flights4 <-
flights3 %>% 
  mutate(dest = fct_lump_prop(dest, prop = .025)) 

flights4 %>% 
  count(dest, sort = T)

flights4 %>% 
  filter(!is_extreme) %>% 
  select(dep_delay, dest, origin, carrier) %>% 
  group_by(dest, origin, carrier) %>% 
  summarise(depdelay_md = median(dep_delay, na.rm = T)) %>% 
  ggplot() +
  aes(x = origin, y = depdelay_md, color = origin) +
  facet_grid(dest ~ carrier) +
  geom_point()
  

flights4 %>% 
  filter(!is_extreme) %>% 
  select(dep_delay, dest, origin, carrier) %>% 
  group_by(dest, origin, carrier) %>% 
  summarise(depdelay_md = median(dep_delay, na.rm = T))
#> # A tibble: 128 x 4
#> # Groups:   dest, origin [37]
#>    dest  origin carrier depdelay_md
#>    <fct> <chr>  <chr>         <dbl>
#>  1 ATL   EWR    9E               -6
#>  2 ATL   EWR    DL               -3
#>  3 ATL   EWR    EV               -2
#>  4 ATL   EWR    UA               -1
#>  5 ATL   JFK    9E               -2
#>  6 ATL   JFK    DL               -1
#>  7 ATL   LGA    DL               -3
#>  8 ATL   LGA    EV               30
#>  9 ATL   LGA    FL                0
#> 10 ATL   LGA    MQ               -4
#> # … with 118 more rows

flights4_sum <- 
  flights4 %>% 
  filter(!is_extreme) %>% 
  select(month, origin, dep_delay) %>% 
  drop_na() %>% 
  group_by(month, origin) %>% 
  summarise(delay_md = median(dep_delay),
            delay_iqr = IQR(dep_delay),
            delay_n = n()) %>% 
  mutate(month = factor(month),
         delay_n = as.numeric(delay_n))
  

flights4 %>% 
  filter(!is_extreme) %>% 
  select(month, origin, dep_delay) %>% 
  mutate(month = factor(month)) %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = month, y = dep_delay, color = origin) +
  geom_violin() +
  geom_point(data = flights4_sum, 
             aes(y = delay_md,
                 x = month)) +
  facet_wrap( ~ origin)

8 Fazit

Keine heiße Spur bisher. Allerdings erlaubt explorative Datenanalyse nur die gleichzeitige Betrachtung von zwei bis drei, vielleicht vier Variablen. Sind die Zusammenhänge komplizierter in dem Sinne, dass mehrere Variablen für einen Effekt zusammenwirken, so ist es mit explorativen Methoden schwer zu finden.

Der nächste logische Schritt: Wir müssen modellieren. Mit Modellierungsmethoden lassen sich auch hochdimensionale Zusammenhänge finden.

9 Achtung

Diese Analyse ist rein explorativ in dem Sinne, dass keine Hypothesen getestet werden. Es ist damit zu rechnen, dass falsch-positive Befunde auftauchen. Alle Ergebnisse sollten anhand neuer Daten validiert werden.

10 Reproduzierbarkeit

#> ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.1.0 (2021-05-18)
#>  os       macOS Big Sur 10.16         
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       Europe/Berlin               
#>  date     2021-06-22                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
#>  package     * version  date       lib source                        
#>  assertthat    0.2.1    2019-03-21 [1] CRAN (R 4.1.0)                
#>  backports     1.2.1    2020-12-09 [1] CRAN (R 4.1.0)                
#>  blogdown      1.3      2021-04-14 [2] CRAN (R 4.1.0)                
#>  bookdown      0.22     2021-04-22 [2] CRAN (R 4.1.0)                
#>  broom         0.7.6    2021-04-05 [1] CRAN (R 4.1.0)                
#>  bslib         0.2.5.1  2021-05-18 [1] CRAN (R 4.1.0)                
#>  cachem        1.0.5    2021-05-15 [2] CRAN (R 4.1.0)                
#>  callr         3.7.0    2021-04-20 [1] CRAN (R 4.1.0)                
#>  cellranger    1.1.0    2016-07-27 [1] CRAN (R 4.1.0)                
#>  cli           2.5.0    2021-04-26 [1] CRAN (R 4.1.0)                
#>  codetools     0.2-18   2020-11-04 [2] CRAN (R 4.1.0)                
#>  colorspace    2.0-1    2021-05-04 [1] CRAN (R 4.1.0)                
#>  crayon        1.4.1    2021-02-08 [1] CRAN (R 4.1.0)                
#>  data.table    1.14.0   2021-02-21 [1] CRAN (R 4.1.0)                
#>  DBI           1.1.1    2021-01-15 [1] CRAN (R 4.1.0)                
#>  dbplyr        2.1.1    2021-04-06 [1] CRAN (R 4.1.0)                
#>  desc          1.3.0    2021-03-05 [2] CRAN (R 4.1.0)                
#>  devtools      2.4.1    2021-05-05 [2] CRAN (R 4.1.0)                
#>  digest        0.6.27   2020-10-24 [1] CRAN (R 4.1.0)                
#>  dplyr       * 1.0.6    2021-05-05 [1] CRAN (R 4.1.0)                
#>  ellipsis      0.3.2    2021-04-29 [1] CRAN (R 4.1.0)                
#>  evaluate      0.14     2019-05-28 [1] CRAN (R 4.1.0)                
#>  fansi         0.5.0    2021-05-25 [1] CRAN (R 4.1.0)                
#>  farver        2.1.0    2021-02-28 [1] CRAN (R 4.1.0)                
#>  fastDummies * 1.6.3    2020-11-29 [1] CRAN (R 4.1.0)                
#>  fastmap       1.1.0    2021-01-25 [2] CRAN (R 4.1.0)                
#>  forcats     * 0.5.1    2021-01-27 [1] CRAN (R 4.1.0)                
#>  fs            1.5.0    2020-07-31 [1] CRAN (R 4.1.0)                
#>  generics      0.1.0    2020-10-31 [1] CRAN (R 4.1.0)                
#>  ggplot2     * 3.3.3    2020-12-30 [1] CRAN (R 4.1.0)                
#>  glue          1.4.2    2020-08-27 [1] CRAN (R 4.1.0)                
#>  gtable        0.3.0    2019-03-25 [1] CRAN (R 4.1.0)                
#>  haven         2.4.1    2021-04-23 [1] CRAN (R 4.1.0)                
#>  highr         0.9      2021-04-16 [1] CRAN (R 4.1.0)                
#>  hms           1.1.0    2021-05-17 [1] CRAN (R 4.1.0)                
#>  htmltools     0.5.1.1  2021-01-22 [1] CRAN (R 4.1.0)                
#>  httr          1.4.2    2020-07-20 [1] CRAN (R 4.1.0)                
#>  jquerylib     0.1.4    2021-04-26 [1] CRAN (R 4.1.0)                
#>  jsonlite      1.7.2    2020-12-09 [1] CRAN (R 4.1.0)                
#>  knitr         1.33     2021-04-24 [1] CRAN (R 4.1.0)                
#>  labeling      0.4.2    2020-10-20 [1] CRAN (R 4.1.0)                
#>  lifecycle     1.0.0    2021-02-15 [1] CRAN (R 4.1.0)                
#>  lubridate     1.7.10   2021-02-26 [1] CRAN (R 4.1.0)                
#>  magrittr      2.0.1    2020-11-17 [1] CRAN (R 4.1.0)                
#>  memoise       2.0.0    2021-01-26 [2] CRAN (R 4.1.0)                
#>  modelr        0.1.8    2020-05-19 [1] CRAN (R 4.1.0)                
#>  munsell       0.5.0    2018-06-12 [1] CRAN (R 4.1.0)                
#>  pillar        1.6.1    2021-05-16 [1] CRAN (R 4.1.0)                
#>  pkgbuild      1.2.0    2020-12-15 [2] CRAN (R 4.1.0)                
#>  pkgconfig     2.0.3    2019-09-22 [1] CRAN (R 4.1.0)                
#>  pkgload       1.2.1    2021-04-06 [2] CRAN (R 4.1.0)                
#>  prettyunits   1.1.1    2020-01-24 [1] CRAN (R 4.1.0)                
#>  printr      * 0.1.1    2021-01-27 [1] CRAN (R 4.1.0)                
#>  processx      3.5.2    2021-04-30 [1] CRAN (R 4.1.0)                
#>  ps            1.6.0    2021-02-28 [1] CRAN (R 4.1.0)                
#>  purrr       * 0.3.4    2020-04-17 [1] CRAN (R 4.1.0)                
#>  R6            2.5.0    2020-10-28 [1] CRAN (R 4.1.0)                
#>  Rcpp          1.0.6    2021-01-15 [1] CRAN (R 4.1.0)                
#>  readr       * 1.4.0    2020-10-05 [1] CRAN (R 4.1.0)                
#>  readxl        1.3.1    2019-03-13 [1] CRAN (R 4.1.0)                
#>  remotes       2.4.0    2021-06-02 [2] CRAN (R 4.1.0)                
#>  reprex        2.0.0    2021-04-02 [1] CRAN (R 4.1.0)                
#>  rlang         0.4.11   2021-04-30 [1] CRAN (R 4.1.0)                
#>  rmarkdown     2.8      2021-05-07 [1] CRAN (R 4.1.0)                
#>  rprojroot     2.0.2    2020-11-15 [2] CRAN (R 4.1.0)                
#>  rstudioapi    0.13     2020-11-12 [1] CRAN (R 4.1.0)                
#>  rvest         1.0.0    2021-03-09 [1] CRAN (R 4.1.0)                
#>  sass          0.4.0    2021-05-12 [1] CRAN (R 4.1.0)                
#>  scales        1.1.1    2020-05-11 [1] CRAN (R 4.1.0)                
#>  sessioninfo   1.1.1    2018-11-05 [2] CRAN (R 4.1.0)                
#>  stringi       1.6.2    2021-05-17 [1] CRAN (R 4.1.0)                
#>  stringr     * 1.4.0    2019-02-10 [1] CRAN (R 4.1.0)                
#>  testthat      3.0.2    2021-02-14 [2] CRAN (R 4.1.0)                
#>  tibble      * 3.1.2    2021-05-16 [1] CRAN (R 4.1.0)                
#>  tidyr       * 1.1.3    2021-03-03 [1] CRAN (R 4.1.0)                
#>  tidyselect    1.1.1    2021-04-30 [1] CRAN (R 4.1.0)                
#>  tidyverse   * 1.3.1    2021-04-15 [1] CRAN (R 4.1.0)                
#>  usethis       2.0.1    2021-02-10 [2] CRAN (R 4.1.0)                
#>  utf8          1.2.1    2021-03-12 [1] CRAN (R 4.1.0)                
#>  vctrs         0.3.8    2021-04-29 [1] CRAN (R 4.1.0)                
#>  withr         2.4.2    2021-04-18 [1] CRAN (R 4.1.0)                
#>  xfun          0.23     2021-05-15 [1] CRAN (R 4.1.0)                
#>  xml2          1.3.2    2020-04-23 [1] CRAN (R 4.1.0)                
#>  yaml          2.2.1.99 2021-06-18 [1] Github (viking/r-yaml@4788abe)
#> 
#> [1] /Users/sebastiansaueruser/Library/R/x86_64/4.1/library
#> [2] /Library/Frameworks/R.framework/Versions/4.1/Resources/library