data > opinion

Tom Alby

4 Explorative Datenanalyse mit dem tidyverse


Sie sind hier: start / lehrveranstaltungen / datenanalyse data science machine learning / 04 eda tidyverse /

Explorative Datenanalysen mit dem tidyverse

Das tidyverse

Das tidyverse ist eine Erweiterung zu R, die Explorative Datenanalysen stark vereinfacht. Das tidyverse stammt von Hadley Wickham, dem Chief Data Scientist von RStudio. Scherzhafterweise wird das tidyverse von Hadleys Anhängern auch manchmal als Hadleyverse bezeichnet.

DasTidyverse enthält mehrere Libraries:

  • ggplot2: Erweiterte Plotting-Möglichkeiten
  • dplyr: Eine beliebte Bibliothek zur Daten-Manipulation, die im Kurs ausführlicher besprochen wird
  • tidyr: Zum “Aufräumen” von Daten
  • readr: Zum Einlesen von Daten
  • purrr: Erweiterung der Funktionalen Programmierung in R
  • tibble: Eine moderne Version der Dataframes; wird das tidyverse geladen, dann werden aus Dataframes automatisch Tibbles
  • stringr: Library zur Bearbeitung von Strings
  • forcats: Library zur Lösung von Problemen mit Factors

Im Kurs werden nicht alle Pakete behandelt. Es ist möglich, dass auch nur einzelne Libraries installiert werden, mit install.packages(tidyverse) werden alle Libraries auf einmal installiert.

#install.packages("tidyverse")
library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.3     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

dplyr

Im Kurs werden wir uns vor allem mit dplyr beschäftigen (das wird Deutsch “dieplaier” ausgesprochen). Eine wunderbare Einführung gibt der Meister himself hier:

Für diejenigen, die nicht das ganze Video sehen wollen, hier die Kernaussagen. dplyr verwendet 5 Haupt-Verben, mit denen grundlegende Befehle ausgeführt werden.

  • select: Zum Auswählen von Variablen
  • filter: Zum Auswählen von Zeilen
  • mutate: Zum Erzeugen von neuen Variablen/Spalten oder auch zum Verändern dieser
  • summarize: Ergebnisse berechnen (im Video sagt Hadley, dass zu einer Zeile reduziert wird, aber man kann auch mehrere Zeilen als Ergebnis haben)
  • arrange: Zum Sortieren von Zeilen

Ein wichtiger Punkt im Tidyverse ist die so genannte Pipe. Unter UNIX sind Pipes ein beliebtes Mittel, um die Ausgabe aus einer Funktion als Input für eine andere Funktion zu nutzen, z.B.

grep ‘2018’ bericht.csv | wc -l

In diesem Beispiel unter UNIX wird in der Datei bericht.csv nach allen Zeilen gesucht, in denen der String 2018 vorkommt.Die Ausgabe wären dann genau diese Zeilen. Mit dem Operator | wird diese Ausgabe in den nächsten Befehl geleitet. wc -l (steht für word count mit dem Parameter lines, so dass Zeilen und nicht Wörter gezählt werden) zählt die Zeilen. Wir erfahren also mit diesem Befehl, wie viele Zeilen den String 2018 enthalten. Pipes sind ein sehr effizientes Mittel unter UNIX. Genau diesen Mechanismus verwendet das tidyverse, wobei hier der Fokus auf der Lesbarkeit liegt.

In dplyr wird anstatt des “|” das Pipe-Symbol %>% verwendet (früher %.% wie oben im Video zu sehen). Wir nutzen für die Demonstration von dplyr wieder den Datensatz ChickWeight. Wieder schauen wir uns dazu das durchschnittliche Gewicht der Hühner an, die mit der Diät Nr. 1 ernährt wurden, dieses Mal verwenden wir aber die dplyr-Syntax:

ChickWeight %>%
  filter(Diet==1) %>%
  filter(Time == 21) %>%
  summary(.)
##      weight           Time        Chick    Diet  
##  Min.   : 96.0   Min.   :21   13     : 1   1:16  
##  1st Qu.:137.5   1st Qu.:21   9      : 1   2: 0  
##  Median :166.0   Median :21   20     : 1   3: 0  
##  Mean   :177.8   Mean   :21   10     : 1   4: 0  
##  3rd Qu.:207.5   3rd Qu.:21   17     : 1         
##  Max.   :305.0   Max.   :21   19     : 1         
##                               (Other):10

Zunächst einmal “werfen” wir den Datensatz auf die nächste Funktion filter, die dann nur die Zeilen auswirft, bei denen die Hühner die Diät Nr. 1 erhalten haben. In der nächsten Zeile wird der Output des vorherigen Filters auf den nächsten Filter geworfen, bei dem wir nur die Datenreihen auswählen, in denen der letzte Zeitpunkt des Wiegens aufgezeichnet wurde. Das Ergebnis “werfen” wir dann auf den Befehl summary(). Damit der Befehl das Ergebnis aus dem vorherigen Schritt empfangen kann, wird der Punkt in der Klammer benötigt. Vergleichen wir das noch mal mit dem Befehl, den wir oben verwendet hatten:

summary(ChickWeight$weight[ChickWeight$Diet == 1 & ChickWeight$Time == 21])

Ja, kann man machen, und mit etwas Erfahrung sieht man auch sehr schnell, was dieser eine Ausdruck bedeutet, aber gerade für Anfänger ist er um einiges komplexer und fehleranfälliger. Auf der anderen Seite kann man so in einer Zeile ausdrücken, wofür man in dplyr vier benötigt. Dennoch ist der Code so schneller nachvollziehbar.

Nun einmal im Vergleich dazu die Hühner, die Diät Nr. 4 bekommen haben:

ChickWeight %>%
  filter(Diet==4) %>%
  filter(Time == 21) %>%
  summary(.)
##      weight           Time        Chick   Diet 
##  Min.   :196.0   Min.   :21   45     :1   1:0  
##  1st Qu.:204.0   1st Qu.:21   43     :1   2:0  
##  Median :237.0   Median :21   41     :1   3:0  
##  Mean   :238.6   Mean   :21   47     :1   4:9  
##  3rd Qu.:264.0   3rd Qu.:21   49     :1        
##  Max.   :322.0   Max.   :21   46     :1        
##                               (Other):3

Offensichtlich hat die Diät Nr. 4 dazu geführt, dass die Hühner schwerer geworden sind. Ob dieser Unterschied statistisch signifikant ist, wird zu einem späteren Zeitpunkt genauer beleuchtet.

Diese wenigen Verben helfen uns also, Datensätze zu filtern, sortieren und auch zu manipulieren. Der besondere Vorteil von dyplr ist hier, dass es eine etwas einfachere Herangehensweise an R bietet.

Stichprobenverteilung des Mittelwerts

Nun wird es einmal kurz etwas komplizierter. Wir gehen dazu ein kleines Gedankenspiel durch: Wir stellen uns vor, dass wir ein Sample beziehungsweise eine Stichprobe nehmen, und von dieser Stichprobe nehmen wir den Durchschnitt.

x <- rnorm(1000000, 3, .25)
mean(sample(x, 50, replace = FALSE, prob = NULL))
## [1] 2.999553

Der Durchschnitt dieses Samples ist wahrscheinlich nicht genau bei dem richtigen Durchschnitt. Und jetzt bedienen wir uns eines Tricks: Wir nehmen nicht ein Sample, sondern ganz ganz viele, in diesem Fall 1.000, und bei jedem Sample schreiben wir uns den Durchschnitt auf. Und dann plotten wir diese Durchschnitte in ein Histogramm:

n <- 1000
sample_means = rep(NA, n)
for(i in 1:n){
  sample_means[i] = mean(sample(x, 50, replace = FALSE, prob = NULL))
}
hist(sample_means, breaks=100, main = "Histogramm von Sample Means", xlab = "Sample Means")

Sieht aus wie eine Normalverteilung, was zunächst überraschend erscheint, aber eigentlich ganz logisch ist. Denn es ist sehr wahrscheinlich, dass wir die Durchschnitte der durch Zufall gezogenen Stichproben eher zu dem tatsächlichen Durchschnitt tendieren als zu einem Durchschnitt weiter draußen.

Das Tolle an diesem kleinen Experiment ist aber, dass wir auch eine Normalverteilung bekommen, wenn wir mehrfach das Mean aus einer Nicht-Normalverteilung ziehen, wie hier zum Beispiel aus einer rechtsschiefen Verteilung. Diese basteln wir uns erst einmal:

e <- rgamma(10000, 3)
hist(e, breaks = 100)

Und nun ziehen wir wieder 1000 Samples daraus:

n <- 1000
sample_means = rep(NA, n)
for(i in 1:n){
  sample_means[i] = mean(sample(e, 50, replace = FALSE, prob = NULL))
}
hist(sample_means, breaks=100, main = "Histogramm von Sample Means", xlab = "Sample Means")

Diese sogenannte Stichprobenverteilung des Mittelwerts klingt zunächst einmal etwas akademisch, aber sie hat einen tollen Zweck. Denn eine Normalverteilung hat die Eigenschaft, wie wir oben gesehen haben, dass 95% der Werte innerhalb von +/- 1.96 Standardabweichungen liegen. Das heißt, wenn wir ein Sample nehmen und den Durchschnitt berechnen, dann liegt dieser ermittelte Wert mit einer Wahrscheinlichkeit von 95% innerhalb von +/- 1.96 Standardabweichungen des tatsächlichen Durchschnitts. Nur 5% der Werte liegen außerhalb des Bereichs.

Wenn ich nun also ein Sample nehme, dann weiß ich eventuell nicht, wie der Durchschnitt der Population ist. Aber ich weiß, egal welchen Durchschnitt ich berechne aus meinem Sample, wie hoch die Wahrscheinlichkeit ist, dass er aus dem Bereich innerhalb von +/- 1.96 Standardabweichungen stammt.

Der t-Test

Zurück zu unseren Hühnern: Wir haben vier Hühnergruppen, die jeweils unterschiedliche Diäten bekommen haben. Wenn wir wissen wollen, ob die Unterschiede in der Zunahme statistisch signifikant sind, dann ist eines der populärsten Verfahren der t-Test. In einer Variante dieses Tests werden zwei Stichproben genommen und deren Mittelwert miteinander verglichen (normalerweise müsste das untersuchte Merkmal normalverteilt sein, was hier nicht der Fall ist, aber das klammern wir jetzt einmal aus). Erfunden wurde dieser Test von William Sealy Gosset, der in der Guiness Brauerei in Dublin arbeitete, aber leider sein Verfahren zur Qualitätssicherung nicht veröffentlichen durfte. Daher wählte er das Pseudonym “Student”, so dass der Test auch manchmal Student’s t-Test genannt wird.

Wir hatten uns oben schon die Durchschnitte des Gewichts der Hühner mit zwei unterschiedlichen Diäten angesehen und festgestellt, dass diese sehr unterschiedlich sind. Aber sind sie so unterschiedlich, dass wir daran glauben, dass der Unterschied statistisch signifikant ist? Was bedeutet das überhaupt? Auch hier bedienen wir uns wieder eines Tricks. Wir sagen nämlich zunächst einmal, dass es keinen Unterschied gibt (Nullhypothese). Wir lehnen diese Nullhypothese aber ab und nehmen die Alternativhypothese (es gibt einen Unterschied an), wenn die Mittelwerte mit einem Konfidenzniveau von 95% nicht aus derselben Population kommen. Anders gesagt: Wir haben zwei Mittelwerte, und wir nehmen die vorher beschriebene Stichprobenverteilung des Mittelwerts zur Hand und schauen dann, ob einer der Werte mindestens 1.96 Standardabweichungen (in diesem Fall nennen wir das Standardfehler) von dem anderen entfernt liegt. Ist das der Fall, dann gehen wir davon aus, dass diese Mittelwerte von Stichproben aus verschiedenen Populationen kommen.

Mit anderen Worten: “Statistisch signifikant” bedeutet nicht, dass die Stichproben definitiv aus unterschiedlichen Populationen kommen. Es könnte zum Beispiel sein, dass wir zwei Stichproben erwischt haben, deren Durchschnitte durch Zufall sehr weit auseinander liegen. Und daher wird immer von einem Konfidenzniveau von 95% gesprochen (die dem 95%-Bereich der Stichprobenverteilung des Mittelwerts entspricht) oder, umgekehrt, weil 100% - 95%, einem Signifikanzniveau von 5%. Diese Grenze kann man sich übrigens selber setzen. Wir hätten auch 99% sagen können. Oder 98%. Es war Ronald Fisher, der die Grenze von 95% vorgeschlagen hatte in seinem Klassiker “Statistical Methods for Research Workers” von 1925:

“The value for which P = .05, or 1 in 20, is 1.96 or nearly 2 ; it is convenient to take this point as a limit in judging whether a deviation is to be considered significant or not.”

Nun rechnen wir das einmal mit 95%. Wir holen uns erst die einzelnen Werte der Diäten 1 und 4:

(diet_vier <- ChickWeight %>%
  filter(Diet==4) %>%
  filter(Time == 21) %>%
  select(weight))
##   weight
## 1    204
## 2    281
## 3    200
## 4    196
## 5    238
## 6    205
## 7    322
## 8    237
## 9    264
(diet_eins <- ChickWeight %>%
  filter(Diet==1) %>%
  filter(Time == 21) %>%
  select(weight))
##    weight
## 1     205
## 2     215
## 3     202
## 4     157
## 5     223
## 6     157
## 7     305
## 8      98
## 9     124
## 10    175
## 11    205
## 12     96
## 13    266
## 14    142
## 15    157
## 16    117

Und nun führen wir den t-Test durch:

t.test(diet_eins,diet_vier)
## 
##  Welch Two Sample t-test
## 
## data:  diet_eins and diet_vier
## t = -2.9525, df = 21.064, p-value = 0.007587
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -103.62720  -17.98391
## sample estimates:
## mean of x mean of y 
##  177.7500  238.5556

Der p-Wert ist nun der Wert, auf den wir achten müssen, denn er besagt, dass wir hier unter dem Signifikanzniveau von 5% liegen (0.007587 bedeutet 0.7587%!). In diesem Fall ist der Unterschied also statistisch signifikant. Keine Sorge übrigens, wenn das mit dem p-Wert nicht so einfach zu verstehen ist, denn selbst Wissenschaftler haben Probleme mit dem p-Wert.

Flüge aus New York und die NAs

Wir verwenden für unsere nächste Analyse wieder den Datensatz nycflights13. Uns interessiert, ob bestimmte Routen, Airlines oder andere Faktoren häufiger Verspätungen haben als andere. Zunächst laden wir die Library:

library(nycflights13)

Dann schauen wir uns die verschiedenen Datensätze an:

flights
## # A tibble: 336,776 x 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # … with 336,766 more rows, and 11 more variables: arr_delay <dbl>,
## #   carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>
summary(flights)
##       year          month             day           dep_time    sched_dep_time
##  Min.   :2013   Min.   : 1.000   Min.   : 1.00   Min.   :   1   Min.   : 106  
##  1st Qu.:2013   1st Qu.: 4.000   1st Qu.: 8.00   1st Qu.: 907   1st Qu.: 906  
##  Median :2013   Median : 7.000   Median :16.00   Median :1401   Median :1359  
##  Mean   :2013   Mean   : 6.549   Mean   :15.71   Mean   :1349   Mean   :1344  
##  3rd Qu.:2013   3rd Qu.:10.000   3rd Qu.:23.00   3rd Qu.:1744   3rd Qu.:1729  
##  Max.   :2013   Max.   :12.000   Max.   :31.00   Max.   :2400   Max.   :2359  
##                                                  NA's   :8255                 
##    dep_delay          arr_time    sched_arr_time   arr_delay       
##  Min.   : -43.00   Min.   :   1   Min.   :   1   Min.   : -86.000  
##  1st Qu.:  -5.00   1st Qu.:1104   1st Qu.:1124   1st Qu.: -17.000  
##  Median :  -2.00   Median :1535   Median :1556   Median :  -5.000  
##  Mean   :  12.64   Mean   :1502   Mean   :1536   Mean   :   6.895  
##  3rd Qu.:  11.00   3rd Qu.:1940   3rd Qu.:1945   3rd Qu.:  14.000  
##  Max.   :1301.00   Max.   :2400   Max.   :2359   Max.   :1272.000  
##  NA's   :8255      NA's   :8713                  NA's   :9430      
##    carrier              flight       tailnum             origin         
##  Length:336776      Min.   :   1   Length:336776      Length:336776     
##  Class :character   1st Qu.: 553   Class :character   Class :character  
##  Mode  :character   Median :1496   Mode  :character   Mode  :character  
##                     Mean   :1972                                        
##                     3rd Qu.:3465                                        
##                     Max.   :8500                                        
##                                                                         
##      dest              air_time        distance         hour      
##  Length:336776      Min.   : 20.0   Min.   :  17   Min.   : 1.00  
##  Class :character   1st Qu.: 82.0   1st Qu.: 502   1st Qu.: 9.00  
##  Mode  :character   Median :129.0   Median : 872   Median :13.00  
##                     Mean   :150.7   Mean   :1040   Mean   :13.18  
##                     3rd Qu.:192.0   3rd Qu.:1389   3rd Qu.:17.00  
##                     Max.   :695.0   Max.   :4983   Max.   :23.00  
##                     NA's   :9430                                  
##      minute        time_hour                  
##  Min.   : 0.00   Min.   :2013-01-01 05:00:00  
##  1st Qu.: 8.00   1st Qu.:2013-04-04 13:00:00  
##  Median :29.00   Median :2013-07-03 10:00:00  
##  Mean   :26.23   Mean   :2013-07-03 05:22:54  
##  3rd Qu.:44.00   3rd Qu.:2013-10-01 07:00:00  
##  Max.   :59.00   Max.   :2013-12-31 23:00:00  
## 
airports
## # A tibble: 1,458 x 8
##    faa   name                       lat    lon   alt    tz dst   tzone          
##    <chr> <chr>                    <dbl>  <dbl> <dbl> <dbl> <chr> <chr>          
##  1 04G   Lansdowne Airport         41.1  -80.6  1044    -5 A     America/New_Yo…
##  2 06A   Moton Field Municipal A…  32.5  -85.7   264    -6 A     America/Chicago
##  3 06C   Schaumburg Regional       42.0  -88.1   801    -6 A     America/Chicago
##  4 06N   Randall Airport           41.4  -74.4   523    -5 A     America/New_Yo…
##  5 09J   Jekyll Island Airport     31.1  -81.4    11    -5 A     America/New_Yo…
##  6 0A9   Elizabethton Municipal …  36.4  -82.2  1593    -5 A     America/New_Yo…
##  7 0G6   Williams County Airport   41.5  -84.5   730    -5 A     America/New_Yo…
##  8 0G7   Finger Lakes Regional A…  42.9  -76.8   492    -5 A     America/New_Yo…
##  9 0P2   Shoestring Aviation Air…  39.8  -76.6  1000    -5 U     America/New_Yo…
## 10 0S9   Jefferson County Intl     48.1 -123.    108    -8 A     America/Los_An…
## # … with 1,448 more rows
airlines
## # A tibble: 16 x 2
##    carrier name                       
##    <chr>   <chr>                      
##  1 9E      Endeavor Air Inc.          
##  2 AA      American Airlines Inc.     
##  3 AS      Alaska Airlines Inc.       
##  4 B6      JetBlue Airways            
##  5 DL      Delta Air Lines Inc.       
##  6 EV      ExpressJet Airlines Inc.   
##  7 F9      Frontier Airlines Inc.     
##  8 FL      AirTran Airways Corporation
##  9 HA      Hawaiian Airlines Inc.     
## 10 MQ      Envoy Air                  
## 11 OO      SkyWest Airlines Inc.      
## 12 UA      United Air Lines Inc.      
## 13 US      US Airways Inc.            
## 14 VX      Virgin America             
## 15 WN      Southwest Airlines Co.     
## 16 YV      Mesa Airlines Inc.
weather
## # A tibble: 26,115 x 15
##    origin  year month   day  hour  temp  dewp humid wind_dir wind_speed
##    <chr>  <int> <int> <int> <int> <dbl> <dbl> <dbl>    <dbl>      <dbl>
##  1 EWR     2013     1     1     1  39.0  26.1  59.4      270      10.4 
##  2 EWR     2013     1     1     2  39.0  27.0  61.6      250       8.06
##  3 EWR     2013     1     1     3  39.0  28.0  64.4      240      11.5 
##  4 EWR     2013     1     1     4  39.9  28.0  62.2      250      12.7 
##  5 EWR     2013     1     1     5  39.0  28.0  64.4      260      12.7 
##  6 EWR     2013     1     1     6  37.9  28.0  67.2      240      11.5 
##  7 EWR     2013     1     1     7  39.0  28.0  64.4      240      15.0 
##  8 EWR     2013     1     1     8  39.9  28.0  62.2      250      10.4 
##  9 EWR     2013     1     1     9  39.9  28.0  62.2      260      15.0 
## 10 EWR     2013     1     1    10  41    28.0  59.6      260      13.8 
## # … with 26,105 more rows, and 5 more variables: wind_gust <dbl>, precip <dbl>,
## #   pressure <dbl>, visib <dbl>, time_hour <dttm>
planes
## # A tibble: 3,322 x 9
##    tailnum  year type          manufacturer   model  engines seats speed engine 
##    <chr>   <int> <chr>         <chr>          <chr>    <int> <int> <int> <chr>  
##  1 N10156   2004 Fixed wing m… EMBRAER        EMB-1…       2    55    NA Turbo-…
##  2 N102UW   1998 Fixed wing m… AIRBUS INDUST… A320-…       2   182    NA Turbo-…
##  3 N103US   1999 Fixed wing m… AIRBUS INDUST… A320-…       2   182    NA Turbo-…
##  4 N104UW   1999 Fixed wing m… AIRBUS INDUST… A320-…       2   182    NA Turbo-…
##  5 N10575   2002 Fixed wing m… EMBRAER        EMB-1…       2    55    NA Turbo-…
##  6 N105UW   1999 Fixed wing m… AIRBUS INDUST… A320-…       2   182    NA Turbo-…
##  7 N107US   1999 Fixed wing m… AIRBUS INDUST… A320-…       2   182    NA Turbo-…
##  8 N108UW   1999 Fixed wing m… AIRBUS INDUST… A320-…       2   182    NA Turbo-…
##  9 N109UW   1999 Fixed wing m… AIRBUS INDUST… A320-…       2   182    NA Turbo-…
## 10 N110UW   1999 Fixed wing m… AIRBUS INDUST… A320-…       2   182    NA Turbo-…
## # … with 3,312 more rows

Das erste, was wir uns anschauen, ist, ob Flugzeit und Ankunftsverspätung miteinander korrelieren:

flights %>%
  select(arr_delay,air_time) %>%
  plot(.)

Auf dem Plot sieht das schon mal nicht so aus, führen wir nun einen Korrelationstest durch:

flights %>%
  select(arr_delay,air_time) %>%
  cor
##           arr_delay air_time
## arr_delay         1       NA
## air_time         NA        1

Hier hat anscheinend etwas nicht funktioniert. Schauen wir uns die Daten einmal genauer an, so sehen wir, dass es darin NAs gibt:

flights %>%
  filter(!is.na(air_time))
## # A tibble: 327,346 x 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # … with 327,336 more rows, and 11 more variables: arr_delay <dbl>,
## #   carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>

Damit kann die Funktion nicht umgehen, daher müssen wir diese Werte exkludieren. Dies geschieht indem wir filtern nach is.na(), diesen Ausdruck aber negieren, indem wir ein ! davor schreiben. Und da auch in der Spalte air_time NAs vorkommen können, tun wir dies dort auch. Da Programmierer faul sind, schreiben wir das hier nicht in zwei Filter-Zeilen, sondern schreiben die beiden Filter-Ausdrücke in eine Zeile:

flights %>%
  select(arr_delay,air_time) %>%
  filter(!is.na(arr_delay) || !is.na(air_time)) %>%
  cor
##           arr_delay air_time
## arr_delay         1       NA
## air_time         NA        1

Wie wir sehen, sehen wir keine Korrelation. Das “|” oder ||" steht für ein “oder”. Man könnte diese Zeile so interpretieren: “Filter alle Zeilen heraus, in denen entweder in der Spalte arr_delay oder der Spalte air_time ein NA drin ist”. Wobei das nicht ganz stimmt. Das soll das folgende Beispiel verdeutlichen:

x <- data.frame("X1" = c(5,NA,8,NA), "X2" = c(6,7,NA,NA))
x %>%
  filter(!is.na(X1) || !is.na(X2))
##   X1 X2
## 1  5  6
## 2 NA  7
## 3  8 NA
## 4 NA NA

Wird “||” genutzt, so wird also nicht wirklich ein “oder” eingesetzt, denn die 4. Zeile dürfte dann nicht ausgespuckt werden. Dies wird deutlich, wenn wir nur ein “|” verwenden:

x %>%
  filter(!is.na(X1) | !is.na(X2))
##   X1 X2
## 1  5  6
## 2 NA  7
## 3  8 NA

Wir haben also zwei Arten von “oder”, und dasselbe gilt für & und &&, die “und” bedeuten. Beim || wird abgebrochen, wenn bereits das erste Element im Vergleich TRUE ist.

flights %>%
  filter(!is.na(arr_delay)) %>%
  filter(!is.na(dep_delay)) %>%
  select(arr_delay,dep_delay) %>%
  cor
##           arr_delay dep_delay
## arr_delay 1.0000000 0.9148028
## dep_delay 0.9148028 1.0000000

Verknüpfen von Daten

Wir haben oben mehrere verschiedene Data Frames, die jeweils in einer Beziehung zueinander stehen. So stehen die richtigen Namen der Reiseziele in dem Data Frame airports, in dem Data Frame stehen nur die Kurzbezeichnungen. Die Kurz-Codes der Airlines im flights Data Frame sind auch nicht unbedingt lesbarer, werden aber im Data Frame airlines übersetzt. Im ersten Schritt sollen genau diese beiden Data Frames miteinander verknüpft werden. Dies geschieht mit dem Befehl left_join. Zusätzlich entfernen wir noch die ursprüngliche Spalte und bennen die zugefügte Spalte um:

flights
## # A tibble: 336,776 x 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # … with 336,766 more rows, and 11 more variables: arr_delay <dbl>,
## #   carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>
airlines
## # A tibble: 16 x 2
##    carrier name                       
##    <chr>   <chr>                      
##  1 9E      Endeavor Air Inc.          
##  2 AA      American Airlines Inc.     
##  3 AS      Alaska Airlines Inc.       
##  4 B6      JetBlue Airways            
##  5 DL      Delta Air Lines Inc.       
##  6 EV      ExpressJet Airlines Inc.   
##  7 F9      Frontier Airlines Inc.     
##  8 FL      AirTran Airways Corporation
##  9 HA      Hawaiian Airlines Inc.     
## 10 MQ      Envoy Air                  
## 11 OO      SkyWest Airlines Inc.      
## 12 UA      United Air Lines Inc.      
## 13 US      US Airways Inc.            
## 14 VX      Virgin America             
## 15 WN      Southwest Airlines Co.     
## 16 YV      Mesa Airlines Inc.
(my_flights <- flights %>%
  left_join(airlines, by = "carrier") %>%
  select(-carrier) %>%
  rename(carrier = name))
## # A tibble: 336,776 x 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # … with 336,766 more rows, and 11 more variables: arr_delay <dbl>,
## #   flight <int>, tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>,
## #   distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>, carrier <chr>

Nicht immer aber haben die Spalten denselben Namen in den Data Frames, die man miteinander verbinden möchte. Hier muss left_join geholfen werden, anhand welchen Keys die Data Frames miteinander verbunden werden sollen.

(my_flights <- my_flights %>%
  left_join(airports, c("dest" = "faa")) %>%
   rename(airport = name))
## # A tibble: 336,776 x 26
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # … with 336,766 more rows, and 18 more variables: arr_delay <dbl>,
## #   flight <int>, tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>,
## #   distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>, carrier <chr>,
## #   airport <chr>, lat <dbl>, lon <dbl>, alt <dbl>, tz <dbl>, dst <chr>,
## #   tzone <chr>

Nachdem das nun erledigt ist, schauen wir uns an, ob es eine Häufung der Verspätungen bei bestimmten Flughäfen gibt:

my_flights %>%
  select(airport,arr_delay) %>%
  filter(!is.na(arr_delay)) %>%
  filter(arr_delay > 0) %>%
  group_by(airport) %>%
  mutate(mean_delay = mean(arr_delay)) %>%
  mutate(median_delay = median(arr_delay)) %>%
  mutate(max_delay = max(arr_delay)) %>%
  mutate(sum_flights = n()) %>%
  select(airport,mean_delay,median_delay,max_delay,sum_flights) %>%
  unique() %>%
  arrange(desc(mean_delay))
## # A tibble: 100 x 5
## # Groups:   airport [100]
##    airport                         mean_delay median_delay max_delay sum_flights
##    <chr>                                <dbl>        <dbl>     <dbl>       <int>
##  1 Cherry Capital Airport                66.7         40         220          35
##  2 Tulsa Intl                            60.4         43         262         193
##  3 Birmingham Intl                       60.1         44         291         121
##  4 Mc Ghee Tyson                         59.9         39         281         307
##  5 Des Moines Intl                       57.0         35.5       322         260
##  6 Will Rogers World                     56.6         42         262         201
##  7 Columbia Metropolitan                 53.6         36         224          87
##  8 Richmond Intl                         53.3         31         463        1196
##  9 Dane Co Rgnl Truax Fld                52.8         34         364         285
## 10 Cincinnati Northern Kentucky I…       52.4         29         989        1680
## # … with 90 more rows

Das gleiche tun wir noch mal mit Airlines:

my_flights %>%
  select(carrier,arr_delay) %>%
  filter(!is.na(arr_delay)) %>%
  filter(arr_delay > 0) %>%
  group_by(carrier) %>%
  mutate(mean_delay = mean(arr_delay)) %>%
  mutate(median_delay = median(arr_delay)) %>%
  mutate(max_delay = max(arr_delay)) %>%
  mutate(sum_flights = n()) %>%
  select(carrier,mean_delay,median_delay,max_delay,sum_flights) %>%
  unique() %>%
  arrange(desc(mean_delay))
## # A tibble: 16 x 5
## # Groups:   carrier [16]
##    carrier                     mean_delay median_delay max_delay sum_flights
##    <chr>                            <dbl>        <dbl>     <dbl>       <int>
##  1 SkyWest Airlines Inc.             60.6         47         157          10
##  2 Mesa Airlines Inc.                51.1         27.5       381         258
##  3 Endeavor Air Inc.                 49.3         27         744        6637
##  4 ExpressJet Airlines Inc.          48.3         28         577       24484
##  5 Frontier Airlines Inc.            47.6         24         834         392
##  6 Virgin America                    43.8         17         676        1746
##  7 AirTran Airways Corporation       41.1         19         572        1895
##  8 Southwest Airlines Co.            40.7         19         453        5304
##  9 JetBlue Airways                   40.0         22         497       23609
## 10 American Airlines Inc.            38.3         19        1007       10706
## 11 Envoy Air                         37.9         20        1127       11693
## 12 Delta Air Lines Inc.              37.7         17         931       16413
## 13 United Air Lines Inc.             36.7         19         455       22222
## 14 Hawaiian Airlines Inc.            35.0         14        1272          97
## 15 Alaska Airlines Inc.              34.4         17         198         189
## 16 US Airways Inc.                   29.0         15         492        7349

Ein großer Nachteil dieser Tabelle ist, dass man das Wesen der Daten nicht schnell erfassen kann. Da wir bei den Airlines eine geringe Anzahl von Merkmalsattributen haben, können wir die Daten auch in einem Boxplot visualisieren:

ggplot(data = my_flights, mapping = aes(x = carrier, y = arr_delay)) +
  theme(axis.text.x=element_text(angle=90,hjust=1)) +
                    geom_boxplot() 
## Warning: Removed 9430 rows containing non-finite values (stat_boxplot).

Nun schauen wir uns noch an, ob das Flugzeugmodell einen Einfluss auf die Verspätung hat. Zunächst einmal verbinden wir die Datensätze:

(my_flights <- my_flights %>%
  left_join(planes, by = "tailnum"))
## # A tibble: 336,776 x 34
##    year.x month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##     <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1   2013     1     1      517            515         2      830            819
##  2   2013     1     1      533            529         4      850            830
##  3   2013     1     1      542            540         2      923            850
##  4   2013     1     1      544            545        -1     1004           1022
##  5   2013     1     1      554            600        -6      812            837
##  6   2013     1     1      554            558        -4      740            728
##  7   2013     1     1      555            600        -5      913            854
##  8   2013     1     1      557            600        -3      709            723
##  9   2013     1     1      557            600        -3      838            846
## 10   2013     1     1      558            600        -2      753            745
## # … with 336,766 more rows, and 26 more variables: arr_delay <dbl>,
## #   flight <int>, tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>,
## #   distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>, carrier <chr>,
## #   airport <chr>, lat <dbl>, lon <dbl>, alt <dbl>, tz <dbl>, dst <chr>,
## #   tzone <chr>, year.y <int>, type <chr>, manufacturer <chr>, model <chr>,
## #   engines <int>, seats <int>, speed <int>, engine <chr>
my_flights %>%
  group_by(manufacturer) %>%
  filter(arr_delay > 0) %>%
  select(manufacturer,arr_delay) %>%
  mutate(mean_delay = mean(arr_delay)) %>%
  mutate(median_delay = median(arr_delay)) %>%
  mutate(max_delay = max(arr_delay)) %>%
  mutate(sum_flights = n()) %>%
  select(manufacturer,mean_delay, median_delay,max_delay,sum_flights) %>%
  unique() %>%
  arrange(desc(mean_delay))
## # A tibble: 36 x 5
## # Groups:   manufacturer [36]
##    manufacturer       mean_delay median_delay max_delay sum_flights
##    <chr>                   <dbl>        <dbl>     <dbl>       <int>
##  1 AVIAT AIRCRAFT INC       86           97         147           5
##  2 AGUSTA SPA               76.1         29         316          17
##  3 FRIEDEMANN JON           63.0         35         279          20
##  4 LAMBERT RICHARD          61           35.5       246          18
##  5 CANADAIR                 52.9         24.5       421         674
##  6 BOMBARDIER INC           50.2         28         744       10614
##  7 LEBLANC GLENN T          48.4         36         163           8
##  8 BELL                     47.5         14.5       307          28
##  9 EMBRAER                  45.6         26         538       29126
## 10 HURLEY JAMES LARRY       44.7         15         165           6
## # … with 26 more rows

Jetzt schauen wir uns an, ob das Wetter einen Einfluss hat:

my_flights %>%
  left_join(weather, by="time_hour")
## # A tibble: 1,006,987 x 48
##    year.x month.x day.x dep_time sched_dep_time dep_delay arr_time
##     <int>   <int> <int>    <int>          <int>     <dbl>    <int>
##  1   2013       1     1      517            515         2      830
##  2   2013       1     1      517            515         2      830
##  3   2013       1     1      517            515         2      830
##  4   2013       1     1      533            529         4      850
##  5   2013       1     1      533            529         4      850
##  6   2013       1     1      533            529         4      850
##  7   2013       1     1      542            540         2      923
##  8   2013       1     1      542            540         2      923
##  9   2013       1     1      542            540         2      923
## 10   2013       1     1      544            545        -1     1004
## # … with 1,006,977 more rows, and 41 more variables: sched_arr_time <int>,
## #   arr_delay <dbl>, flight <int>, tailnum <chr>, origin.x <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour.x <dbl>, minute <dbl>,
## #   time_hour <dttm>, carrier <chr>, airport <chr>, lat <dbl>, lon <dbl>,
## #   alt <dbl>, tz <dbl>, dst <chr>, tzone <chr>, year.y <int>, type <chr>,
## #   manufacturer <chr>, model <chr>, engines <int>, seats <int>, speed <int>,
## #   engine <chr>, origin.y <chr>, year <int>, month.y <int>, day.y <int>,
## #   hour.y <int>, temp <dbl>, dewp <dbl>, humid <dbl>, wind_dir <dbl>,
## #   wind_speed <dbl>, wind_gust <dbl>, precip <dbl>, pressure <dbl>,
## #   visib <dbl>

Was ist das Problem mit diesem Join? Wir haben plötzlich über 1 Million Datensätze, vorher hatten wir etwas mehr als 336.000. Ein Problem ist, dass das verbindende Glied nicht allein die Zeit sein kann, denn wir haben 3 Flughäfen, und da diese weit auseinander liegen, kann das Wetter von einem zum nächsten Flughafen sehr unterschiedlich ausfallen. left_join kommt aber auch damit klar, wenn man keinen Parameter angibt; dann werden die Spalten ausgesucht, die matchen:

(my_flights <- my_flights %>%
  left_join(weather))
## Joining, by = c("month", "day", "origin", "hour", "time_hour")
## # A tibble: 336,776 x 44
##    year.x month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##     <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1   2013     1     1      517            515         2      830            819
##  2   2013     1     1      533            529         4      850            830
##  3   2013     1     1      542            540         2      923            850
##  4   2013     1     1      544            545        -1     1004           1022
##  5   2013     1     1      554            600        -6      812            837
##  6   2013     1     1      554            558        -4      740            728
##  7   2013     1     1      555            600        -5      913            854
##  8   2013     1     1      557            600        -3      709            723
##  9   2013     1     1      557            600        -3      838            846
## 10   2013     1     1      558            600        -2      753            745
## # … with 336,766 more rows, and 36 more variables: arr_delay <dbl>,
## #   flight <int>, tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>,
## #   distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>, carrier <chr>,
## #   airport <chr>, lat <dbl>, lon <dbl>, alt <dbl>, tz <dbl>, dst <chr>,
## #   tzone <chr>, year.y <int>, type <chr>, manufacturer <chr>, model <chr>,
## #   engines <int>, seats <int>, speed <int>, engine <chr>, year <int>,
## #   temp <dbl>, dewp <dbl>, humid <dbl>, wind_dir <dbl>, wind_speed <dbl>,
## #   wind_gust <dbl>, precip <dbl>, pressure <dbl>, visib <dbl>

Die Formel für die Umrechnung von Fahrenheit in Celsius sieht wie folgt aus: C = (°F − 32) / 1,8. Man könnte sie wie folgt umrechnen:

(my_flights <- my_flights %>%
  mutate(celsius = (temp-32) / 1.8))
## # A tibble: 336,776 x 45
##    year.x month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##     <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1   2013     1     1      517            515         2      830            819
##  2   2013     1     1      533            529         4      850            830
##  3   2013     1     1      542            540         2      923            850
##  4   2013     1     1      544            545        -1     1004           1022
##  5   2013     1     1      554            600        -6      812            837
##  6   2013     1     1      554            558        -4      740            728
##  7   2013     1     1      555            600        -5      913            854
##  8   2013     1     1      557            600        -3      709            723
##  9   2013     1     1      557            600        -3      838            846
## 10   2013     1     1      558            600        -2      753            745
## # … with 336,766 more rows, and 37 more variables: arr_delay <dbl>,
## #   flight <int>, tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>,
## #   distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>, carrier <chr>,
## #   airport <chr>, lat <dbl>, lon <dbl>, alt <dbl>, tz <dbl>, dst <chr>,
## #   tzone <chr>, year.y <int>, type <chr>, manufacturer <chr>, model <chr>,
## #   engines <int>, seats <int>, speed <int>, engine <chr>, year <int>,
## #   temp <dbl>, dewp <dbl>, humid <dbl>, wind_dir <dbl>, wind_speed <dbl>,
## #   wind_gust <dbl>, precip <dbl>, pressure <dbl>, visib <dbl>, celsius <dbl>

Wir plotten nun verschiedene Faktoren und schauen, welche davon Einfluss auf Verspätungen haben. Dieser Plot kann ein wenig dauern :)

my_flights %>%
  filter(arr_delay > 30) %>%
  #filter(!is.na(arr_delay) || !is.na(celsius) || !is.na(visib))
  select(arr_delay,celsius,visib,precip,wind_speed,wind_gust,dep_delay) %>%
  plot(.)

Der größte Faktor für Verspätungen scheint die verspätete Abflugszeit zu sein, Überraschung :) Das Wetter scheint ansonsten keinen Einfluss zu haben.

Da wir nun einen transformierten Datensatz haben und diese Schritte nicht jedes Mal durchgehen wollen, speichern wir unseren neuen Datensatz zum Schluß einmal ab:

write_csv(my_flights,path="my_flights.csv")

Lineare Regression: Was beeinflusst die Anzahl von täglichen Flügen?

Das folgende Beispiel stammt vom Meister himself und seinem Co-Autor, as dem Buch R for Data Science. Hierzu werden zwei weitere Libraries benötigt.

library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:dplyr':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(modelr)

Es wird derselbe Datensatz verwendet, in diesem Fall aber schauen wir uns nur an, wie viele Flüge es pro Wochentag gegeben hat:

daily <- flights %>% 
  mutate(date = make_date(year, month, day)) %>% 
  group_by(date) %>% 
  summarise(n = n())
daily
## # A tibble: 365 x 2
##    date           n
##    <date>     <int>
##  1 2013-01-01   842
##  2 2013-01-02   943
##  3 2013-01-03   914
##  4 2013-01-04   915
##  5 2013-01-05   720
##  6 2013-01-06   832
##  7 2013-01-07   933
##  8 2013-01-08   899
##  9 2013-01-09   902
## 10 2013-01-10   932
## # … with 355 more rows

Dies plotten wir nun mit ggplot2:

ggplot(daily, aes(date, n)) + 
  geom_line()

Offensichtlich sehen wir hier nicht so viel, außer dass es anscheinend jede Woche einen Peak nach unten gibt.

daily <- daily %>% 
  mutate(wday = wday(date, label = TRUE))
ggplot(daily, aes(wday, n)) + 
  geom_boxplot()

Am Wochenende wird weniger geflogen, vor allem Samstags, wie Wickham feststellt. Man fliegt vielleicht Sonntag Abends eher, weil man Montag morgens irgendwo einen Termin hat, aber der Samstag ist im Vergleich zu den anderen Tagen ruhiger.

Wir würden aber gerne verstehen, ob es noch andere Muster sieht, die unter diesem Muster des Samstags verborgen sind. Dazu bilden wir erst einmal ein Model, wir nutzen dazu eine lineare Regression. Dabei haben wir zwei Variablen, die erklärende Variable (der Wochentag) und die abhängige Variable (die Anzahl der Flüge). Eine lineare Regression funktioniert vereinfacht so, dass der Algorithmus versucht, eine Gerade durch die Punkte so zu legen, dass die geringsten quadrierten Abweichungen entstehen. Wieso quadriert? Wir hatten diesen Ansatz schon bei der Standardabweichung gesehen, nämlich dass wir quadrieren, um das negative Vorzeichen loszuwerden. Denn eine Abweichung kann auch hier negativ sein. Diese Methode der kleinsten Quadrate wird uns in ähnlicher Form wieder begegnen.

mod <- lm(n ~ wday, data = daily)

grid <- daily %>% 
  data_grid(wday) %>% 
  add_predictions(mod, "n")

ggplot(daily, aes(wday, n)) + 
  geom_boxplot() +
  geom_point(data = grid, colour = "red", size = 4)

Im nächsten Schritt werden die Residuen berechnet und dann visualisiert. Residuen sind die Abweichungen von dem Wert, den das Modell errechnet hat. Diese plotten wir dann:

daily <- daily %>% 
  add_residuals(mod)
daily %>% 
  ggplot(aes(date, resid)) + 
  geom_ref_line(h = 0) + 
  geom_line()

Wir sehen nun, an welchen Tagen wir unerwartete Ausreißer haben. Dies können wir auch noch mal in einer Tabelle darstellen:

daily %>% 
  filter(resid < -100)
## # A tibble: 11 x 4
##    date           n wday  resid
##    <date>     <int> <ord> <dbl>
##  1 2013-01-01   842 Tue   -109.
##  2 2013-01-20   786 Sun   -105.
##  3 2013-05-26   729 Sun   -162.
##  4 2013-07-04   737 Thu   -229.
##  5 2013-07-05   822 Fri   -145.
##  6 2013-09-01   718 Sun   -173.
##  7 2013-11-28   634 Thu   -332.
##  8 2013-11-29   661 Fri   -306.
##  9 2013-12-24   761 Tue   -190.
## 10 2013-12-25   719 Wed   -244.
## 11 2013-12-31   776 Tue   -175.

Die Gründe für diese Ausreißer sind zum Teil offensichtlich, zum Teil weniger (Wickham hat die weniger offensichtlichen Ausreißer zur Hausaufgabe deklariert :)

Weitere Informationen zum tidyverse

Ein weiteres, informatives Video mit Hadley Wickham ist hier zu finden: