Introduction

Visualisation is an important tool for insight generation, but it is rare that you get the data in exactly the right form you need. Often you’ll need to create some new variables or summaries, or maybe you just want to rename the variables or reorder the observations in order to make the data a little easier to work with. You’ll learn how to do all that (and more!) in this chapter, which will teach you how to transform your data using the dplyr package and a new dataset on flights departing New York City in 2013.

Prerequisites

In this chapter we’re going to focus on how to use the dplyr package, another core member of the tidyverse. We’ll illustrate the key ideas using data from the nycflights13 package, and use ggplot2 to help us understand the data.

Libraries ***

library(nycflights13)
library(tidyverse)
library(dplyr)
library(ggplot2)

Take careful note of the conflicts message that’s printed when you load the tidyverse. It tells you that dplyr overwrites some functions in base R. If you want to use the base version of these functions after loading dplyr, you’ll need to use their full names: stats::filter() and stats::lag().

nycflights13

To explore the basic data manipulation verbs of dplyr, we’ll use nycflights13::flights. This data frame contains all 336,776 flights that departed from New York City in 2013. The data comes from the US Bureau of Transportation Statistics, and is documented in ?flights.

flights
#> # A tibble: 336,776 x 19
#>    year month   day 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      533            529         4      850
#> 3  2013     1     1      542            540         2      923
#> 4  2013     1     1      544            545        -1     1004
#> 5  2013     1     1      554            600        -6      812
#> 6  2013     1     1      554            558        -4      740
#> # ... with 3.368e+05 more rows, and 12 more variables:
#> #   sched_arr_time <int>, 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>

You might notice that this data frame prints a little differently from other data frames you might have used in the past: it only shows the first few rows and all the columns that fit on one screen. (To see the whole dataset, you can run View(flights) which will open the dataset in the RStudio viewer). It prints differently because it’s a tibble. Tibbles are data frames, but slightly tweaked to work better in the tidyverse. For now, you don’t need to worry about the differences; we’ll come back to tibbles in more detail in wrangle.

You might also have noticed the row of three (or four) letter abbreviations under the column names. These describe the type of each variable:

  • int stands for integers.

  • dbl stands for doubles, or real numbers.

  • chr stands for character vectors, or strings.

  • dttm stands for date-times (a date + a time).

There are three other common types of variables that aren’t used in this dataset but you’ll encounter later in the book:

  • lgl stands for logical, vectors that contain only TRUE or FALSE.

  • fctr stands for factors, which R uses to represent categorical variables with fixed possible values.

  • date stands for dates.

dplyr basics

In this chapter you are going to learn the five key dplyr functions that allow you to solve the vast majority of your data manipulation challenges:

  • Pick observations by their values (filter()).
  • Reorder the rows (arrange()).
  • Pick variables by their names (select()).
  • Create new variables with functions of existing variables (mutate()).
  • Collapse many values down to a single summary (summarise()).

These can all be used in conjunction with group_by() which changes the scope of each function from operating on the entire dataset to operating on it group-by-group. These six functions provide the verbs for a language of data manipulation.

All verbs work similarly:

  1. The first argument is a data frame.

  2. The subsequent arguments describe what to do with the data frame, using the variable names (without quotes).

  3. The result is a new data frame.

Together these properties make it easy to chain together multiple simple steps to achieve a complex result. Let’s dive in and see how these verbs work.

Filter rows with filter()

filter() allows you to subset observations based on their values. The first argument is the name of the data frame. The second and subsequent arguments are the expressions that filter the data frame. For example, we can select all flights on January 1st with:

filter(flights, month == 1, day == 1)
#> # A tibble: 842 x 19
#>    year month   day 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      533            529         4      850
#> 3  2013     1     1      542            540         2      923
#> 4  2013     1     1      544            545        -1     1004
#> 5  2013     1     1      554            600        -6      812
#> 6  2013     1     1      554            558        -4      740
#> # ... with 836 more rows, and 12 more variables: sched_arr_time <int>,
#> #   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>

When you run that line of code, dplyr executes the filtering operation and returns a new data frame. dplyr functions never modify their inputs, so if you want to save the result, you’ll need to use the assignment operator, <-:

jan1 <- filter(flights, month == 1, day == 1)

R either prints out the results, or saves them to a variable. If you want to do both, you can wrap the assignment in parentheses:

(dec25 <- filter(flights, month == 12, day == 25))
#> # A tibble: 719 x 19
#>    year month   day dep_time sched_dep_time dep_delay arr_time
#>   <int> <int> <int>    <int>          <int>     <dbl>    <int>
#> 1  2013    12    25      456            500        -4      649
#> 2  2013    12    25      524            515         9      805
#> 3  2013    12    25      542            540         2      832
#> 4  2013    12    25      546            550        -4     1022
#> 5  2013    12    25      556            600        -4      730
#> 6  2013    12    25      557            600        -3      743
#> # ... with 713 more rows, and 12 more variables: sched_arr_time <int>,
#> #   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>

Comparisons

To use filtering effectively, you have to know how to select the observations that you want using the comparison operators. R provides the standard suite: >, >=, <, <=, != (not equal), and == (equal).

When you’re starting out with R, the easiest mistake to make is to use = instead of == when testing for equality. When this happens you’ll get an informative error:

filter(flights, month = 1)
#> Error: `month` (`month = 1`) must not be named, do you need `==`?

There’s another common problem you might encounter when using ==: floating point numbers. These results might surprise you!

sqrt(2) ^ 2 == 2
#> [1] FALSE
1/49 * 49 == 1
#> [1] FALSE
near(), sqrt()

Computers use finite precision arithmetic (they obviously can’t store an infinite number of digits!) so remember that every number you see is an approximation. Instead of relying on ==, use near():

near(sqrt(2) ^ 2,  2)
#> [1] TRUE
near(1 / 49 * 49, 1)
#> [1] TRUE

Logical operators

Multiple arguments to filter() are combined with “and”: every expression must be true in order for a row to be included in the output. For other types of combinations, you’ll need to use Boolean operators yourself: & is “and”, | is “or”, and ! is “not”. Figure @ref(fig:bool-ops) shows the complete set of Boolean operations.

Complete set of boolean operations. `x` is the left-hand circle, `y` is the right-hand circle, and the shaded region show which parts each operator selects.

Complete set of boolean operations. x is the left-hand circle, y is the right-hand circle, and the shaded region show which parts each operator selects.

The following code finds all flights that departed in November or December:

filter(flights, month == 11 | month == 12)

The order of operations doesn’t work like English. You can’t write filter(flights, month == 11 | 12), which you might literally translate into “finds all flights that departed in November or December”. Instead it finds all months that equal 11 | 12, an expression that evaluates to TRUE. In a numeric context (like here), TRUE becomes one, so this finds all flights in January, not November or December. This is quite confusing!

A useful short-hand for this problem is x %in% y. This will select every row where x is one of the values in y. We could use it to rewrite the code above:

nov_dec <- filter(flights, month %in% c(11, 12))

Sometimes you can simplify complicated subsetting by remembering De Morgan’s law: !(x & y) is the same as !x | !y, and !(x | y) is the same as !x & !y. For example, if you wanted to find flights that weren’t delayed (on arrival or departure) by more than two hours, you could use either of the following two filters:

filter(flights, !(arr_delay > 120 | dep_delay > 120))
filter(flights, arr_delay <= 120, dep_delay <= 120)

As well as & and |, R also has && and ||. Don’t use them here! You’ll learn when you should use them in [conditional execution].

Whenever you start using complicated, multipart expressions in filter(), consider making them explicit variables instead. That makes it much easier to check your work. You’ll learn how to create new variables shortly.

Missing values

One important feature of R that can make comparison tricky are missing values, or NAs (“not availables”). NA represents an unknown value so missing values are “contagious”: almost any operation involving an unknown value will also be unknown.

NA > 5
#> [1] NA
10 == NA
#> [1] NA
NA + 10
#> [1] NA
NA / 2
#> [1] NA

The most confusing result is this one:

NA == NA
#> [1] NA

It’s easiest to understand why this is true with a bit more context:

# Let x be Mary's age. We don't know how old she is.
x <- NA

# Let y be John's age. We don't know how old he is.
y <- NA

# Are John and Mary the same age?
x == y
#> [1] NA
# We don't know!
is.na()

If you want to determine if a value is missing, use is.na():

is.na(x)
#> [1] TRUE

filter() only includes rows where the condition is TRUE; it excludes both FALSE and NA values. If you want to preserve missing values, ask for them explicitly:

df <- tibble(x = c(1, NA, 3))
filter(df, x > 1)
#> # A tibble: 1 x 1
#>       x
#>   <dbl>
#> 1     3
filter(df, is.na(x) | x > 1)
#> # A tibble: 2 x 1
#>       x
#>   <dbl>
#> 1    NA
#> 2     3

Exercises

  1. Find all flights that

    1. Had an arrival delay of two or more hours

      filter(flights, arr_delay > 120)
    2. Flew to Houston (IAH or HOU)

      filter(flights, dest == "IAH" | dest == "HOU")
    3. Were operated by United, American, or Delta

      filter(flights, carrier =="AA" | carrier == "UA" | carrier =="DL")
      
      # OR
      
      filter(flights, carrier %in% c("AA", "UA", "DL"))
    4. Departed in summer (July, August, and September)

      filter(flights, month %in% c(7, 8, 9))
    5. Arrived more than two hours late, but didn’t leave late

      filter(flights, arr_delay > 120 & dep_delay <= 0)
    6. Were delayed by at least an hour, but made up over 30 minutes in flight

      filter(flights, dep_delay > 60 & (arr_time < (sched_arr_time + 30)))
    7. Departed between midnight and 6am (inclusive)

      filter(flights, dep_time >= 0000 & dep_time <= 0600)
between()
  1. Another useful dplyr filtering helper is between(). What does it do? Can you use it to simplify the code needed to answer the previous challenges?

    filter(flights, between(dep_time, 0000, 0600))
  2. How many flights have a missing dep_time? What other variables are missing? What might these rows represent?

    filter(flights,is.na(dep_time)) # 8255 flights
    
    # What other variables are missing? dep_delay, arr_time, arr_delay
    
    # What might these rows represent? Cancelled flights
  3. Why is NA ^ 0 not missing? Why is NA | TRUE not missing? Why is FALSE & NA not missing? Can you figure out the general rule? (NA * 0 is a tricky counterexample!)

    NA ^ 0
    # [1] 1
    
    NA | TRUE
    #  [1] TRUE
    
    NA & FALSE
    # [1] FALSE
    
    NA * 0
    # [1] NA

Arrange rows with arrange()

arrange() works similarly to filter() except that instead of selecting rows, it changes their order. It takes a data frame and a set of column names (or more complicated expressions) to order by. If you provide more than one column name, each additional column will be used to break ties in the values of preceding columns:

arrange(flights, year, month, day)
#> # A tibble: 336,776 x 19
#>    year month   day 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      533            529         4      850
#> 3  2013     1     1      542            540         2      923
#> 4  2013     1     1      544            545        -1     1004
#> 5  2013     1     1      554            600        -6      812
#> 6  2013     1     1      554            558        -4      740
#> # ... with 3.368e+05 more rows, and 12 more variables:
#> #   sched_arr_time <int>, 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>

Use desc() to re-order by a column in descending order:

arrange(flights, desc(arr_delay))
#> # A tibble: 336,776 x 19
#>    year month   day dep_time sched_dep_time dep_delay arr_time
#>   <int> <int> <int>    <int>          <int>     <dbl>    <int>
#> 1  2013     1     9      641            900      1301     1242
#> 2  2013     6    15     1432           1935      1137     1607
#> 3  2013     1    10     1121           1635      1126     1239
#> 4  2013     9    20     1139           1845      1014     1457
#> 5  2013     7    22      845           1600      1005     1044
#> 6  2013     4    10     1100           1900       960     1342
#> # ... with 3.368e+05 more rows, and 12 more variables:
#> #   sched_arr_time <int>, 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>

Missing values are always sorted at the end:

df <- tibble(x = c(5, 2, NA))
arrange(df, x)
#> # A tibble: 3 x 1
#>       x
#>   <dbl>
#> 1     2
#> 2     5
#> 3    NA
arrange(df, desc(x))
#> # A tibble: 3 x 1
#>       x
#>   <dbl>
#> 1     5
#> 2     2
#> 3    NA

Exercises

  1. How could you use arrange() to sort all missing values to the start? (Hint: use is.na()).

    arrange(df, desc(is.na(x)))
    
  2. Sort flights to find the most delayed flights. Find the flights that left earliest.

    arrange(flights, desc(dep_delay))
    arrange(flights, dep_delay)
    
  3. Sort flights to find the fastest flights.

    arrange(flights, air_time)
    select(arrange(flights, air_time), flight, month, day, year, air_time)
    
  4. Which flights traveled the longest? Which traveled the shortest?

    select(arrange(flights, desc(air_time)), flight, month, day, year, air_time) 
    select(arrange(flights, air_time), flight, month, day, year, air_time)
    

Select columns with select()

It’s not uncommon to get datasets with hundreds or even thousands of variables. In this case, the first challenge is often narrowing in on the variables you’re actually interested in. select() allows you to rapidly zoom in on a useful subset using operations based on the names of the variables.

select() is not terribly useful with the flights data because we only have 19 variables, but you can still get the general idea:

# Select columns by name
select(flights, year, month, day)
#> # A tibble: 336,776 x 3
#>    year month   day
#>   <int> <int> <int>
#> 1  2013     1     1
#> 2  2013     1     1
#> 3  2013     1     1
#> 4  2013     1     1
#> 5  2013     1     1
#> 6  2013     1     1
#> # ... with 3.368e+05 more rows
# Select all columns between year and day (inclusive)
select(flights, year:day)
#> # A tibble: 336,776 x 3
#>    year month   day
#>   <int> <int> <int>
#> 1  2013     1     1
#> 2  2013     1     1
#> 3  2013     1     1
#> 4  2013     1     1
#> 5  2013     1     1
#> 6  2013     1     1
#> # ... with 3.368e+05 more rows
# Select all columns except those from year to day (inclusive)
select(flights, -(year:day))
#> # A tibble: 336,776 x 16
#>   dep_time sched_dep_time dep_delay arr_time sched_arr_time arr_delay
#>      <int>          <int>     <dbl>    <int>          <int>     <dbl>
#> 1      517            515         2      830            819        11
#> 2      533            529         4      850            830        20
#> 3      542            540         2      923            850        33
#> 4      544            545        -1     1004           1022       -18
#> 5      554            600        -6      812            837       -25
#> 6      554            558        -4      740            728        12
#> # ... with 3.368e+05 more rows, and 10 more variables: carrier <chr>,
#> #   flight <int>, tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>,
#> #   distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>

starts_with(), ends_with(), contains(), matches(), num_range(), one_of()

There are a number of helper functions you can use within select():

  • starts_with("abc"): matches names that begin with “abc”.

  • ends_with("xyz"): matches names that end with “xyz”.

  • contains("ijk"): matches names that contain “ijk”.

  • matches("(.)\\1"): selects variables that match a regular expression. This one matches any variables that contain repeated characters. You’ll learn more about regular expressions in [strings].

  • num_range("x", 1:3) matches x1, x2 and x3.

See ?select for more details.

rename()

select() can be used to rename variables, but it’s rarely useful because it drops all of the variables not explicitly mentioned. Instead, use rename(), which is a variant of select() that keeps all the variables that aren’t explicitly mentioned:

rename(flights, tail_num = tailnum) # new = old
#> # A tibble: 336,776 x 19
#>    year month   day 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      533            529         4      850
#> 3  2013     1     1      542            540         2      923
#> 4  2013     1     1      544            545        -1     1004
#> 5  2013     1     1      554            600        -6      812
#> 6  2013     1     1      554            558        -4      740
#> # ... with 3.368e+05 more rows, and 12 more variables:
#> #   sched_arr_time <int>, arr_delay <dbl>, carrier <chr>, flight <int>,
#> #   tail_num <chr>, origin <chr>, dest <chr>, air_time <dbl>,
#> #   distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>
everything()

Another option is to use select() in conjunction with the everything() helper. This is useful if you have a handful of variables you’d like to move to the start of the data frame.

select(flights, time_hour, air_time, everything())
#> # A tibble: 336,776 x 19
#>             time_hour air_time  year month   day dep_time sched_dep_time
#>                <dttm>    <dbl> <int> <int> <int>    <int>          <int>
#> 1 2013-01-01 05:00:00      227  2013     1     1      517            515
#> 2 2013-01-01 05:00:00      227  2013     1     1      533            529
#> 3 2013-01-01 05:00:00      160  2013     1     1      542            540
#> 4 2013-01-01 05:00:00      183  2013     1     1      544            545
#> 5 2013-01-01 06:00:00      116  2013     1     1      554            600
#> 6 2013-01-01 05:00:00      150  2013     1     1      554            558
#> # ... with 3.368e+05 more rows, and 12 more variables: dep_delay <dbl>,
#> #   arr_time <int>, sched_arr_time <int>, arr_delay <dbl>, carrier <chr>,
#> #   flight <int>, tailnum <chr>, origin <chr>, dest <chr>, distance <dbl>,
#> #   hour <dbl>, minute <dbl>

Exercises

  1. Brainstorm as many ways as possible to select dep_time, dep_delay, arr_time, and arr_delay from flights.

    select(flights, dep_time, dep_delay, arr_time, arr_delay) # YES x4
    select(flights, starts_with("dep"), starts_with("arr")) # YES x4
    select(flights, ends_with("_time"), ends_with("_delay")) #NO x7
    select(flights, contains("arr"), contains("dep")) #NO x7
  2. What happens if you include the name of a variable multiple times in a select() call?

    It is returned once, not mulitple times

  3. What does the one_of() function do? Why might it be helpful in conjunction with this vector?

    vars <- c("year", "month", "day", "dep_delay", "arr_delay")
    vars
    #> [1] "year"      "month"     "day"       "dep_delay" "arr_delay"
    These functions allow you to select variables based on their names:
    • starts_with(): starts with a prefix
    • ends_with(): ends with a prefix
    • contains(): contains a literal string
    • matches(): matches a regular expression
    • num_range(): a numerical range like x01, x02, x03.
    • one_of(): variables in character vector.
    • everything(): all variables.
    select(flights, one_of(vars)) #returns each column in vars
    #> # A tibble: 336,776 x 5
    #>    year month   day dep_delay arr_delay
    #>   <int> <int> <int>     <dbl>     <dbl>
    #> 1  2013     1     1         2        11
    #> 2  2013     1     1         4        20
    #> 3  2013     1     1         2        33
    #> 4  2013     1     1        -1       -18
    #> 5  2013     1     1        -6       -25
    #> 6  2013     1     1        -4        12
    #> # ... with 3.368e+05 more rows
  4. Does the result of running the following code surprise you? How do the select helpers deal with case by default? How can you change that default?

    select(flights, contains("TIME"))
    Yes, it is a surprise that helpers are not case sensitive. How do the select helpers deal with case by default?
    • ignore.case If TRUE, the default, ignores case when matching names. How can you change that default?
    select(flights, contains("TIME", ignore.case = FALSE))
    

Add new variables with mutate()

Besides selecting sets of existing columns, it’s often useful to add new columns that are functions of existing columns. That’s the job of mutate().

mutate() always adds new columns at the end of your dataset so we’ll start by creating a narrower dataset so we can see the new variables. Remember that when you’re in RStudio, the easiest way to see all the columns is View().

flight_sml ***
flights_sml <- select(flights, 
  year:day, 
  ends_with("delay"), 
  distance, 
  air_time
)
mutate(flights_sml,
  gain = arr_delay - dep_delay,
  speed = distance / air_time * 60
)
#> # A tibble: 336,776 x 9
#>    year month   day dep_delay arr_delay distance air_time  gain    speed
#>   <int> <int> <int>     <dbl>     <dbl>    <dbl>    <dbl> <dbl>    <dbl>
#> 1  2013     1     1         2        11     1400      227     9 370.0441
#> 2  2013     1     1         4        20     1416      227    16 374.2731
#> 3  2013     1     1         2        33     1089      160    31 408.3750
#> 4  2013     1     1        -1       -18     1576      183   -17 516.7213
#> 5  2013     1     1        -6       -25      762      116   -19 394.1379
#> 6  2013     1     1        -4        12      719      150    16 287.6000
#> # ... with 3.368e+05 more rows

Note that you can refer to columns that you’ve just created:

mutate(flights_sml,
  gain = arr_delay - dep_delay,
  hours = air_time / 60,
  gain_per_hour = gain / hours
)
#> # A tibble: 336,776 x 10
#>    year month   day dep_delay arr_delay distance air_time  gain    hours
#>   <int> <int> <int>     <dbl>     <dbl>    <dbl>    <dbl> <dbl>    <dbl>
#> 1  2013     1     1         2        11     1400      227     9 3.783333
#> 2  2013     1     1         4        20     1416      227    16 3.783333
#> 3  2013     1     1         2        33     1089      160    31 2.666667
#> 4  2013     1     1        -1       -18     1576      183   -17 3.050000
#> 5  2013     1     1        -6       -25      762      116   -19 1.933333
#> 6  2013     1     1        -4        12      719      150    16 2.500000
#> # ... with 3.368e+05 more rows, and 1 more variables: gain_per_hour <dbl>
transmute()

If you only want to keep the new variables, use transmute():

transmute(flights,
  gain = arr_delay - dep_delay,
  hours = air_time / 60,
  gain_per_hour = gain / hours
)
#> # A tibble: 336,776 x 3
#>    gain    hours gain_per_hour
#>   <dbl>    <dbl>         <dbl>
#> 1     9 3.783333      2.378855
#> 2    16 3.783333      4.229075
#> 3    31 2.666667     11.625000
#> 4   -17 3.050000     -5.573770
#> 5   -19 1.933333     -9.827586
#> 6    16 2.500000      6.400000
#> # ... with 3.368e+05 more rows

Useful creation functions

There are many functions for creating new variables that you can use with mutate(). The key property is that the function must be vectorised: it must take a vector of values as input, return a vector with the same number of values as output. There’s no way to list every possible function that you might use, but here’s a selection of functions that are frequently useful:

Arithmetic operators
  • Arithmetic operators: +, -, *, /, ^. These are all vectorised, using the so called “recycling rules”. If one parameter is shorter than the other, it will be automatically extended to be the same length. This is most useful when one of the arguments is a single number: air_time / 60, hours * 60 + minute, etc.

    Arithmetic operators are also useful in conjunction with the aggregate functions you’ll learn about later. For example, x / sum(x) calculates the proportion of a total, and y - mean(y) computes the difference from the mean.

Modular arithmetic
  • Modular arithmetic: %/% (integer division) and %% (remainder), where x == y * (x %/% y) + (x %% y). Modular arithmetic is a handy tool because it allows you to break integers up into pieces. For example, in the flights dataset, you can compute hour and minute from dep_time with:

    transmute(flights,
      dep_time,
      hour = dep_time %/% 100,
      minute = dep_time %% 100
    )
    #> # A tibble: 336,776 x 3
    #>   dep_time  hour minute
    #>      <int> <dbl>  <dbl>
    #> 1      517     5     17
    #> 2      533     5     33
    #> 3      542     5     42
    #> 4      544     5     44
    #> 5      554     5     54
    #> 6      554     5     54
    #> # ... with 3.368e+05 more rows
Logs: log(), log2(), log10()
  • Logs: log(), log2(), log10(). Logarithms are an incredibly useful transformation for dealing with data that ranges across multiple orders of magnitude. They also convert multiplicative relationships to additive, a feature we’ll come back to in modelling.

    All else being equal, I recommend using log2() because it’s easy to interpret: a difference of 1 on the log scale corresponds to doubling on the original scale and a difference of -1 corresponds to halving.

Offsets: lead() and lag()
  • Offsets: lead() and lag() allow you to refer to leading or lagging values. This allows you to compute running differences (e.g. x - lag(x)) or find when values change (x != lag(x)). They are most useful in conjunction with group_by(), which you’ll learn about shortly.

    (x <- 1:10)
    #>  [1]  1  2  3  4  5  6  7  8  9 10
    lag(x)
    #>  [1] NA  1  2  3  4  5  6  7  8  9
    lead(x)
    #>  [1]  2  3  4  5  6  7  8  9 10 NA
Cumulative and rolling aggregates: cumsum(), cumprod(), cummin(), cummax(), cummean()
  • Cumulative and rolling aggregates: R provides functions for running sums, products, mins and maxes: cumsum(), cumprod(), cummin(), cummax(); and dplyr provides cummean() for cumulative means. If you need rolling aggregates (i.e. a sum computed over a rolling window), try the RcppRoll package.
x ***
```text
x
#>  [1]  1  2  3  4  5  6  7  8  9 10
```

```text
cumsum(x)
#>  [1]  1  3  6 10 15 21 28 36 45 55
```

```text
cummean(x)
#>  [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5
```
  • Logical comparisons, <, <=, >, >=, !=, which you learned about earlier. If you’re doing a complex sequence of logical operations it’s often a good idea to store the interim values in new variables so you can check that each step is working as expected.

  • Ranking: there are a number of ranking functions, but you should start with min_rank(). It does the most usual type of ranking (e.g. 1st, 2nd, 2nd, 4th). The default gives smallest values the small ranks; use desc(x) to give the largest values the smallest ranks.

    y <- c(1, 2, 2, NA, 3, 4)
    min_rank(y)
    #> [1]  1  2  2 NA  4  5
    min_rank(desc(y))
    #> [1]  5  3  3 NA  2  1

    If min_rank() doesn’t do what you need, look at the variants row_number(), dense_rank(), percent_rank(), cume_dist(), ntile(). See their help pages for more details.

    row_number(y)
    #> [1]  1  2  3 NA  4  5
    dense_rank(y)
    #> [1]  1  2  2 NA  3  4
    percent_rank(y)
    #> [1] 0.00 0.25 0.25   NA 0.75 1.00
    cume_dist(y)
    #> [1] 0.2 0.6 0.6  NA 0.8 1.0

Exercises

  1. Currently dep_time and sched_dep_time are convenient to look at, but hard to compute with because they’re not really continuous numbers. Convert them to a more convenient representation of number of minutes since midnight.

    transmute(flights,
              dep_time,
              dep_hour = dep_time %/% 100,
              dep_min = dep_time %% 100,
              dep_mins_since_midnight = (dep_hour * 60) + dep_min,
              sched_dep_time,
              sched_dep_mins_since_midnight = ((sched_dep_time %/% 100) * 60) + sched_dep_time %% 100
              )
  2. Compare air_time with arr_time - dep_time. What do you expect to see? What do you see? What do you need to do to fix it?

    transmute(flights,
          air_time,
          enroute = arr_time - dep_time
          )
    #> # A tibble: 336,776 x 2
    #>   air_time enroute
    #>      <dbl>   <int>
    #> 1      227     313
    #> 2      227     317
    #> 3      160     381
    #> 4      183     460
    #> 5      116     258
    #> 6      150     186
    #> # ... with 3.368e+05 more rows

    To fix it, the arr_time and dep_time must be converted to minutes from midnight. However, they are also given in local time zones, so that must be accounted for. As well, some flights could have departed before midnight and landed after midnight, which must also be accounted for.

    Naive results show air_time most often varies from arr_time - dep_time:

    transmute(flights,
      dep_time = ((dep_time %/% 100) * 60) + (dep_time %% 100),
      arr_time = ((arr_time %/% 100) * 60) + (arr_time %% 100),
      enroute = arr_time - dep_time %% (60*24),
      air_time
      )
    #> # A tibble: 336,776 x 4
    #>   dep_time arr_time enroute air_time
    #>      <dbl>    <dbl>   <dbl>    <dbl>
    #> 1      317      510     193      227
    #> 2      333      530     197      227
    #> 3      342      563     221      160
    #> 4      344      604     260      183
    #> 5      354      492     138      116
    #> 6      354      460     106      150
    #> # ... with 3.368e+05 more rows
  3. Compare dep_time, sched_dep_time, and dep_delay. How would you expect those three numbers to be related?

    
    # One might expect dep_time - sched_dep_time == dep_delay
    
    transmute(flights, dep_time, sched_dep_time, dep_delay, x = dep_time - sched_dep_time, y = x == dep_delay)
    
    # However, since arr_time and dep_time are in clock times rather then minutes, an adjustment would need be made (as above) to effect the correct calculation.
    
  4. Find the 10 most delayed flights using a ranking function. How do you want to handle ties? Carefully read the documentation for min_rank().

    ?min_rank
    # equivalent to rank(ties.method = "min")
    
    x <- as.data.frame(select(flights, dep_delay))
    head(min_rank(desc(x)), 10)   
    #>  [1] 114150 103893 114150 144947 258934 209494 234113 185276 185276 163760
  5. What does 1:3 + 1:10 return? Why?

    1:3 + 1:10
    #> Warning in 1:3 + 1:10: longer object length is not a multiple of shorter
    #> object length
    #>  [1]  2  4  6  5  7  9  8 10 12 11
    # returns a warning message because longer object length is not a multiple of shorter object length
  6. What trigonometric functions does R provide?

    Base R provides the following, with angles in radians cos(x), sin(x), tan(x), acos(x), asin(x), atan(x), atan2(y, x), cospi(x), sinpi(x), tanpi(x)

Grouped summaries with summarise()

The last key verb is summarise(). It collapses a data frame to a single row:

summarise(flights, delay = mean(dep_delay, na.rm = TRUE))
#> # A tibble: 1 x 1
#>      delay
#>      <dbl>
#> 1 12.63907

(We’ll come back to what that na.rm = TRUE means very shortly.)

group_by()

summarise() is not terribly useful unless we pair it with group_by(). This changes the unit of analysis from the complete dataset to individual groups. Then, when you use the dplyr verbs on a grouped data frame they’ll be automatically applied “by group”. For example, if we applied exactly the same code to a data frame grouped by date, we get the average delay per date:

by_day <- group_by(flights, year, month, day)
summarise(by_day, delay = mean(dep_delay, na.rm = TRUE))
#> # A tibble: 365 x 4
#> # Groups:   year, month [?]
#>    year month   day     delay
#>   <int> <int> <int>     <dbl>
#> 1  2013     1     1 11.548926
#> 2  2013     1     2 13.858824
#> 3  2013     1     3 10.987832
#> 4  2013     1     4  8.951595
#> 5  2013     1     5  5.732218
#> 6  2013     1     6  7.148014
#> # ... with 359 more rows

Note: Grouping doesn’t change how the data looks (apart from listing how it’s grouped). It changes how it acts with the other dplyr verbs.

Together group_by() and summarise() provide one of the tools that you’ll use most commonly when working with dplyr: grouped summaries. But before we go any further with this, we need to introduce a powerful new idea: the pipe.

Combine multiple operations with the Pipe

Imagine that we want to explore the relationship between the distance and average delay for each location. Using what you know about dplyr, you might write code like this:

by_dest <- group_by(flights, dest)
delay <- summarise(by_dest,
  count = n(),
  dist = mean(distance, na.rm = TRUE),
  delay = mean(arr_delay, na.rm = TRUE)
)
delay <- filter(delay, count > 20, dest != "HNL")

# It looks like delays increase with distance up to ~750 miles 
# and then decrease. Maybe as flights get longer there's more 
# ability to make up delays in the air?
ggplot(data = delay, mapping = aes(x = dist, y = delay)) +
  geom_point(aes(size = count), alpha = 1/3) +
  geom_smooth(se = FALSE)
#> `geom_smooth()` using method = 'loess'

There are three steps to prepare this data:

  1. Group flights by destination.

  2. Summarize to compute distance, average delay, and number of flights.

  3. Filter to remove noisy points and Honolulu airport, which is almost twice as far away as the next closest airport.

This code is a little frustrating to write because we have to give each intermediate data frame a name, even though we don’t care about it. Naming things is hard, so this slows down our analysis.

There’s another way to tackle the same problem with the pipe, %>%:

delays <- flights %>% 
  group_by(dest) %>% 
  summarise(
    count = n(),
    dist = mean(distance, na.rm = TRUE),
    delay = mean(arr_delay, na.rm = TRUE)
  ) %>% 
  filter(count > 20, dest != "HNL")

This focuses on the transformations, not what’s being transformed, which makes the code easier to read. You can read it as a series of imperative statements: group, then summarize, then filter. As suggested by this reading, a good way to pronounce %>% when reading code is “then”.

Behind the scenes, x %>% f(y) turns into f(x, y), and x %>% f(y) %>% g(z) turns into g(f(x, y), z) and so on. You can use the pipe to rewrite multiple operations in a way that you can read left-to-right, top-to-bottom. We’ll use piping frequently from now on because it considerably improves the readability of code, and we’ll come back to it in more detail in [pipes].

Working with the pipe is one of the key criteria for belonging to the tidyverse. The only exception is ggplot2: it was written before the pipe was discovered. Unfortunately, the next iteration of ggplot2, ggvis, which does use the pipe, isn’t quite ready for prime time yet.

Missing values

You may have wondered about the na.rm argument we used above. What happens if we don’t set it?

flights %>% 
  group_by(year, month, day) %>% 
  summarise(mean = mean(dep_delay))
#> # A tibble: 365 x 4
#> # Groups:   year, month [?]
#>    year month   day  mean
#>   <int> <int> <int> <dbl>
#> 1  2013     1     1    NA
#> 2  2013     1     2    NA
#> 3  2013     1     3    NA
#> 4  2013     1     4    NA
#> 5  2013     1     5    NA
#> 6  2013     1     6    NA
#> # ... with 359 more rows

We get a lot of missing values! That’s because aggregation functions obey the usual rule of missing values: if there’s any missing value in the input, the output will be a missing value. Fortunately, all aggregation functions have an na.rm argument which removes the missing values prior to computation:

flights %>% 
  group_by(year, month, day) %>% 
  summarise(mean = mean(dep_delay, na.rm = TRUE))
#> # A tibble: 365 x 4
#> # Groups:   year, month [?]
#>    year month   day      mean
#>   <int> <int> <int>     <dbl>
#> 1  2013     1     1 11.548926
#> 2  2013     1     2 13.858824
#> 3  2013     1     3 10.987832
#> 4  2013     1     4  8.951595
#> 5  2013     1     5  5.732218
#> 6  2013     1     6  7.148014
#> # ... with 359 more rows

In this case, where missing values represent cancelled flights, we could also tackle the problem by first removing the cancelled flights. We’ll save this dataset so we can reuse in the next few examples.

not_cancelled ***
not_cancelled <- flights %>% 
  filter(!is.na(dep_delay), !is.na(arr_delay))
not_cancelled %>% 
  group_by(year, month, day) %>% 
  summarise(mean = mean(dep_delay))
#> # A tibble: 365 x 4
#> # Groups:   year, month [?]
#>    year month   day      mean
#>   <int> <int> <int>     <dbl>
#> 1  2013     1     1 11.435620
#> 2  2013     1     2 13.677802
#> 3  2013     1     3 10.907778
#> 4  2013     1     4  8.965859
#> 5  2013     1     5  5.732218
#> 6  2013     1     6  7.145959
#> # ... with 359 more rows

Counts

n(), sum(!is.na(x))

Whenever you do any aggregation, it’s always a good idea to include either a count (n()), or a count of non-missing values (sum(!is.na(x))). That way you can check that you’re not drawing conclusions based on very small amounts of data. For example, let’s look at the planes (identified by their tail number) that have the highest average delays:

delays <- not_cancelled %>% 
  group_by(tailnum) %>% 
  summarise(
    delay = mean(arr_delay)
  )

ggplot(data = delays, mapping = aes(x = delay)) + 
  geom_freqpoly(binwidth = 10)

Wow, there are some planes that have an average delay of 5 hours (300 minutes)!

The story is actually a little more nuanced. We can get more insight if we draw a scatterplot of number of flights vs. average delay:

delays <- not_cancelled %>% 
  group_by(tailnum) %>% 
  summarise(
    delay = mean(arr_delay, na.rm = TRUE),
    n = n()
  )

ggplot(data = delays, mapping = aes(x = n, y = delay)) + 
  geom_point(alpha = 1/10)

Not surprisingly, there is much greater variation in the average delay when there are few flights. The shape of this plot is very characteristic: whenever you plot a mean (or other summary) vs. group size, you’ll see that the variation decreases as the sample size increases.

When looking at this sort of plot, it’s often useful to filter out the groups with the smallest numbers of observations, so you can see more of the pattern and less of the extreme variation in the smallest groups. This is what the following code does, as well as showing you a handy pattern for integrating ggplot2 into dplyr flows. It’s a bit painful that you have to switch from %>% to +, but once you get the hang of it, it’s quite convenient.

delays %>% 
  filter(n > 25) %>% 
  ggplot(mapping = aes(x = n, y = delay)) + 
    geom_point(alpha = 1/10)


RStudio tip: a useful keyboard shortcut is Cmd/Ctrl + Shift + P. This resends the previously sent chunk from the editor to the console. This is very convenient when you’re (e.g.) exploring the value of n in the example above. You send the whole block once with Cmd/Ctrl + Enter, then you modify the value of n and press Cmd/Ctrl + Shift + P to resend the complete block.


There’s another common variation of this type of pattern. Let’s look at how the average performance of batters in baseball is related to the number of times they’re at bat. Here I use data from the Lahman package to compute the batting average (number of hits / number of attempts) of every major league baseball player.

When I plot the skill of the batter (measured by the batting average, ba) against the number of opportunities to hit the ball (measured by at bat, ab), you see two patterns:

  1. As above, the variation in our aggregate decreases as we get more data points.

  2. There’s a positive correlation between skill (ba) and opportunities to hit the ball (ab). This is because teams control who gets to play, and obviously they’ll pick their best players.

# Convert to a tibble so it prints nicely
batting <- as_tibble(Lahman::Batting)

batters <- batting %>% 
  group_by(playerID) %>% 
  summarise(
    ba = sum(H, na.rm = TRUE) / sum(AB, na.rm = TRUE),
    ab = sum(AB, na.rm = TRUE)
  )

batters %>% 
  filter(ab > 100) %>% 
  ggplot(mapping = aes(x = ab, y = ba)) +
    geom_point() + 
    geom_smooth(se = FALSE)
#> `geom_smooth()` using method = 'gam'

This also has important implications for ranking. If you naively sort on desc(ba), the people with the best batting averages are clearly lucky, not skilled:

batters %>% 
  arrange(desc(ba))
#> # A tibble: 18,915 x 3
#>    playerID    ba    ab
#>       <chr> <dbl> <int>
#> 1 abramge01     1     1
#> 2 banisje01     1     1
#> 3 bartocl01     1     1
#> 4  bassdo01     1     1
#> 5 berrijo01     1     1
#> 6 birasst01     1     2
#> # ... with 1.891e+04 more rows

You can find a good explanation of this problem at http://varianceexplained.org/r/empirical_bayes_baseball/ and http://www.evanmiller.org/how-not-to-sort-by-average-rating.html.

Useful summary functions

Just using means, counts, and sum can get you a long way, but R provides many other useful summary functions:

mean(), median()
  • Measures of location: we’ve used mean(x), but median(x) is also useful. The mean is the sum divided by the length; the median is a value where 50% of x is above it, and 50% is below it.

    It’s sometimes useful to combine aggregation with logical subsetting. We haven’t talked about this sort of subsetting yet, but you’ll learn more about it in [subsetting].

    not_cancelled %>% 
      group_by(year, month, day) %>% 
      summarise(
        avg_delay1 = mean(arr_delay),
        avg_delay2 = mean(arr_delay[arr_delay > 0]) # the average positive delay
      )
    #> # A tibble: 365 x 5
    #> # Groups:   year, month [?]
    #>    year month   day avg_delay1 avg_delay2
    #>   <int> <int> <int>      <dbl>      <dbl>
    #> 1  2013     1     1  12.651023   32.48156
    #> 2  2013     1     2  12.692888   32.02991
    #> 3  2013     1     3   5.733333   27.66087
    #> 4  2013     1     4  -1.932819   28.30976
    #> 5  2013     1     5  -1.525802   22.55882
    #> 6  2013     1     6   4.236429   24.37270
    #> # ... with 359 more rows
sd(), IQR(), mad()
  • Measures of spread: sd(x), IQR(x), mad(x). The mean squared deviation, or standard deviation or sd for short, is the standard measure of spread. The interquartile range IQR() and median absolute deviation mad(x) are robust equivalents that may be more useful if you have outliers.

    # Why is distance to some destinations more variable than to others?
    not_cancelled %>% 
      group_by(dest) %>% 
      summarise(distance_sd = sd(distance)) %>% 
      arrange(desc(distance_sd))
    #> # A tibble: 104 x 2
    #>    dest distance_sd
    #>   <chr>       <dbl>
    #> 1   EGE   10.542765
    #> 2   SAN   10.350094
    #> 3   SFO   10.216017
    #> 4   HNL   10.004197
    #> 5   SEA    9.977993
    #> 6   LAS    9.907786
    #> # ... with 98 more rows
min(), quantile(), max()
  • Measures of rank: min(x), quantile(x, 0.25), max(x). Quantiles are a generalisation of the median. For example, quantile(x, 0.25) will find a value of x that is greater than 25% of the values, and less than the remaining 75%.

    # When do the first and last flights leave each day?
    not_cancelled %>% 
      group_by(year, month, day) %>% 
      summarise(
        first = min(dep_time),
        last = max(dep_time)
      )
    #> # A tibble: 365 x 5
    #> # Groups:   year, month [?]
    #>    year month   day first  last
    #>   <int> <int> <int> <dbl> <dbl>
    #> 1  2013     1     1   517  2356
    #> 2  2013     1     2    42  2354
    #> 3  2013     1     3    32  2349
    #> 4  2013     1     4    25  2358
    #> 5  2013     1     5    14  2357
    #> 6  2013     1     6    16  2355
    #> # ... with 359 more rows
first(), nth(), last()
  • Measures of position: first(x), nth(x, 2), last(x). These work similarly to x[1], x[2], and x[length(x)] but let you set a default value if that position does not exist (i.e. you’re trying to get the 3rd element from a group that only has two elements). For example, we can find the first and last departure for each day:

    not_cancelled %>% 
      group_by(year, month, day) %>% 
      summarise(
        first_dep = first(dep_time), 
        last_dep = last(dep_time)
      )
    #> # A tibble: 365 x 5
    #> # Groups:   year, month [?]
    #>    year month   day first_dep last_dep
    #>   <int> <int> <int>     <int>    <int>
    #> 1  2013     1     1       517     2356
    #> 2  2013     1     2        42     2354
    #> 3  2013     1     3        32     2349
    #> 4  2013     1     4        25     2358
    #> 5  2013     1     5        14     2357
    #> 6  2013     1     6        16     2355
    #> # ... with 359 more rows

    These functions are complementary to filtering on ranks. Filtering gives you all variables, with each observation in a separate row:

    not_cancelled %>% 
      group_by(year, month, day) %>% 
      mutate(r = min_rank(desc(dep_time))) %>% 
      filter(r %in% range(r))
    #> # A tibble: 770 x 20
    #> # Groups:   year, month, day [365]
    #>    year month   day 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     2356           2359        -3      425
    #> 3  2013     1     2       42           2359        43      518
    #> 4  2013     1     2     2354           2359        -5      413
    #> 5  2013     1     3       32           2359        33      504
    #> 6  2013     1     3     2349           2359       -10      434
    #> # ... with 764 more rows, and 13 more variables: sched_arr_time <int>,
    #> #   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>, r <int>
n(), sum(!is.na()), n_distinct(), count()
  • Counts: You’ve seen n(), which takes no arguments, and returns the size of the current group. To count the number of non-missing values, use sum(!is.na(x)). To count the number of distinct (unique) values, use n_distinct(x).

    # Which destinations have the most carriers?
    not_cancelled %>% 
      group_by(dest) %>% 
      summarise(carriers = n_distinct(carrier)) %>% 
      arrange(desc(carriers))
    #> # A tibble: 104 x 2
    #>    dest carriers
    #>   <chr>    <int>
    #> 1   ATL        7
    #> 2   BOS        7
    #> 3   CLT        7
    #> 4   ORD        7
    #> 5   TPA        7
    #> 6   AUS        6
    #> # ... with 98 more rows

    Counts are so useful that dplyr provides a simple helper if all you want is a count:

    not_cancelled %>% 
      count(dest)
    #> # A tibble: 104 x 2
    #>    dest     n
    #>   <chr> <int>
    #> 1   ABQ   254
    #> 2   ACK   264
    #> 3   ALB   418
    #> 4   ANC     8
    #> 5   ATL 16837
    #> 6   AUS  2411
    #> # ... with 98 more rows

    You can optionally provide a weight variable. For example, you could use this to “count” (sum) the total number of miles a plane flew:

    not_cancelled %>% 
      count(tailnum, wt = distance)
    #> # A tibble: 4,037 x 2
    #>   tailnum      n
    #>     <chr>  <dbl>
    #> 1  D942DN   3418
    #> 2  N0EGMQ 239143
    #> 3  N10156 109664
    #> 4  N102UW  25722
    #> 5  N103US  24619
    #> 6  N104UW  24616
    #> # ... with 4,031 more rows
  • Counts and proportions of logical values: sum(x > 10), mean(y == 0). When used with numeric functions, TRUE is converted to 1 and FALSE to 0. This makes sum() and mean() very useful: sum(x) gives the number of TRUEs in x, and mean(x) gives the proportion.

    # How many flights left before 5am? (these usually indicate delayed
    # flights from the previous day)
    not_cancelled %>% 
      group_by(year, month, day) %>% 
      summarise(n_early = sum(dep_time < 500))
    #> # A tibble: 365 x 4
    #> # Groups:   year, month [?]
    #>    year month   day n_early
    #>   <int> <int> <int>   <int>
    #> 1  2013     1     1       0
    #> 2  2013     1     2       3
    #> 3  2013     1     3       4
    #> 4  2013     1     4       3
    #> 5  2013     1     5       3
    #> 6  2013     1     6       2
    #> # ... with 359 more rows
    
    # What proportion of flights are delayed by more than an hour?
    not_cancelled %>% 
      group_by(year, month, day) %>% 
      summarise(hour_perc = mean(arr_delay > 60))
    #> # A tibble: 365 x 4
    #> # Groups:   year, month [?]
    #>    year month   day  hour_perc
    #>   <int> <int> <int>      <dbl>
    #> 1  2013     1     1 0.07220217
    #> 2  2013     1     2 0.08512931
    #> 3  2013     1     3 0.05666667
    #> 4  2013     1     4 0.03964758
    #> 5  2013     1     5 0.03486750
    #> 6  2013     1     6 0.04704463
    #> # ... with 359 more rows

Grouping by multiple variables

When you group by multiple variables, each summary peels off one level of the grouping. That makes it easy to progressively roll up a dataset:

daily ***
daily <- group_by(flights, year, month, day)
(per_day   <- summarise(daily, flights = n()))
#> # A tibble: 365 x 4
#> # Groups:   year, month [?]
#>    year month   day flights
#>   <int> <int> <int>   <int>
#> 1  2013     1     1     842
#> 2  2013     1     2     943
#> 3  2013     1     3     914
#> 4  2013     1     4     915
#> 5  2013     1     5     720
#> 6  2013     1     6     832
#> # ... with 359 more rows
(per_month <- summarise(per_day, flights = sum(flights)))
#> # A tibble: 12 x 3
#> # Groups:   year [?]
#>    year month flights
#>   <int> <int>   <int>
#> 1  2013     1   27004
#> 2  2013     2   24951
#> 3  2013     3   28834
#> 4  2013     4   28330
#> 5  2013     5   28796
#> 6  2013     6   28243
#> # ... with 6 more rows
(per_year  <- summarise(per_month, flights = sum(flights)))
#> # A tibble: 1 x 2
#>    year flights
#>   <int>   <int>
#> 1  2013  336776

Be careful when progressively rolling up summaries: it’s OK for sums and counts, but you need to think about weighting means and variances, and it’s not possible to do it exactly for rank-based statistics like the median. In other words, the sum of groupwise sums is the overall sum, but the median of groupwise medians is not the overall median.

Ungrouping

If you need to remove grouping, and return to operations on ungrouped data, use ungroup().

daily %>% 
  ungroup() %>%             # no longer grouped by date
  summarise(flights = n())  # all flights
#> # A tibble: 1 x 1
#>   flights
#>     <int>
#> 1  336776

Exercises

  1. Brainstorm at least 5 different ways to assess the typical delay characteristics of a group of flights. Consider the following scenarios:

    • A flight is 15 minutes early 50% of the time, and 15 minutes late 50% of the time.
    # First, let's start with a simple sanity check. Do the results look logical when we look at a summary of number of flights for a given flight(number) and carrier comination as compared to how many of them respectivelty are 15 minutes early or late. 
    # Yes, the results look logical. 
    
    not_cancelled %>%
      group_by(flight, carrier) %>%
      summarize(
        n = n(),
        n_early_15 = sum(arr_delay <= -15),
        n_late_15  = sum(arr_delay >= 15)
        )
    #> # A tibble: 5,706 x 5
    #> # Groups:   flight [?]
    #>   flight carrier     n n_early_15 n_late_15
    #>    <int>   <chr> <int>      <int>     <int>
    #> 1      1      AA   361        166        47
    #> 2      1      B6   273         60        88
    #> 3      1      DL    30         18         2
    #> 4      1      UA     3          2         0
    #> 5      1      WN    30          1        10
    #> 6      2      B6    50          4         6
    #> # ... with 5,700 more rows
    # That works, so now let's look at the results in terms of proportions:
    
    # Limit to at least 20 flights to remove small samples.
    # Group by carrier, flight and destination to eliminate overlapping designations.
    # Flights which arrived more than 15 minutes or more early 50% of the time, OR 15 minutes or more late 50% of the time.
    # (No flights averaged both.)
    
    not_cancelled %>%
      group_by(carrier, flight, dest)%>%
      summarize(
        n = n(),
        pct_early_15 = mean(arr_delay <= -15),
        pct_late_15  = mean(arr_delay >= 15)
      ) %>%
      filter(n >= 20, pct_early_15 >= .5 | pct_late_15 >= .5)
    #> # A tibble: 458 x 6
    #> # Groups:   carrier, flight [450]
    #>   carrier flight  dest     n pct_early_15 pct_late_15
    #>     <chr>  <int> <chr> <int>        <dbl>       <dbl>
    #> 1      9E   2907   BUF    44    0.5000000   0.2500000
    #> 2      9E   2918   IND    58    0.6206897   0.2068966
    #> 3      9E   2934   MCI    28    0.6428571   0.1785714
    #> 4      9E   2939   PIT    28    0.5000000   0.1071429
    #> 5      9E   3300   DCA    40    0.6250000   0.1500000
    #> 6      9E   3302   MSP    26    0.5769231   0.1153846
    #> # ... with 452 more rows
    • A flight is always 10 minutes late.
    not_cancelled %>%
      group_by(carrier, flight, dest)%>%
      summarize(
        n = n(),
        pct_late_10 = mean(arr_delay >= 10)
      ) %>%
      filter(n >= 5, pct_late_10 == 1)
    #> # A tibble: 10 x 5
    #> # Groups:   carrier, flight [10]
    #>   carrier flight  dest     n pct_late_10
    #>     <chr>  <int> <chr> <int>       <dbl>
    #> 1      9E   3585   PHL     6           1
    #> 2      EV   3827   CAE     5           1
    #> 3      EV   5058   ATL     7           1
    #> 4      EV   5525   MHT     5           1
    #> 5      EV   5679   CLT     6           1
    #> 6      UA    374   ORD     6           1
    #> # ... with 4 more rows
    • A flight is 30 minutes early 50% of the time, and 30 minutes late 50% of the time.
    not_cancelled %>%
      group_by(carrier, flight, dest)%>%
      summarize(
        n = n(),
        pct_early_30 = mean(arr_delay <= -30),
        pct_late_30 = mean(arr_delay >= 30)
      ) %>%
      filter(n >= 5, pct_early_30 >= .3, pct_late_30 >= .3)
    #> # A tibble: 9 x 6
    #> # Groups:   carrier, flight [9]
    #>   carrier flight  dest     n pct_early_30 pct_late_30
    #>     <chr>  <int> <chr> <int>        <dbl>       <dbl>
    #> 1      9E   4191   BWI    13    0.3076923   0.3076923
    #> 2      EV   4159   CLE     8    0.3750000   0.3750000
    #> 3      UA    271   LAS     6    0.5000000   0.3333333
    #> 4      UA    810   PBI     5    0.4000000   0.4000000
    #> 5      UA   1252   SFO     6    0.5000000   0.3333333
    #> 6      UA   1426   IAH     6    0.3333333   0.3333333
    #> # ... with 3 more rows
    • 99% of the time a flight is on time. 1% of the time it’s 2 hours late.
    not_cancelled %>%
      group_by(carrier, flight, dest) %>%
      summarize(
        n = n(),
        on_time   = mean(arr_delay <= 0),
        late_2hrs = mean(arr_delay >= 120)
      ) %>%
      filter(n >20, on_time >= .99 | late_2hrs >= .01)
    #> # A tibble: 1,838 x 6
    #> # Groups:   carrier, flight [1,729]
    #>   carrier flight  dest     n   on_time  late_2hrs
    #>     <chr>  <int> <chr> <int>     <dbl>      <dbl>
    #> 1      9E   2901   BOS    55 0.6181818 0.05454545
    #> 2      9E   2902   BOS    53 0.6981132 0.01886792
    #> 3      9E   2903   BOS    27 0.6296296 0.03703704
    #> 4      9E   2904   BOS    52 0.5576923 0.03846154
    #> 5      9E   2906   BUF    40 0.5750000 0.07500000
    #> 6      9E   2907   BUF    44 0.7045455 0.02272727
    #> # ... with 1,832 more rows
    # Explore/check:
    flight_2901_9E <- not_cancelled %>%
      filter(carrier =='9E', flight == 2901) %>%
      arrange(desc(arr_delay))
    
    # Yes, the top three rows show that 3 out of 55 flights were more than 2 hours late = 5.45%
    
    flight_2901_9E
    #> # A tibble: 55 x 19
    #>    year month   day dep_time sched_dep_time dep_delay arr_time
    #>   <int> <int> <int>    <int>          <int>     <dbl>    <int>
    #> 1  2013    12    14     1425            825       360     1604
    #> 2  2013    12     5     1126            825       181     1240
    #> 3  2013    12    10     1049            825       144     1208
    #> 4  2013    12    26     1026            830       116     1135
    #> 5  2013    12    18      841            825        16     1054
    #> 6  2013    12    15      914            825        49     1049
    #> # ... with 49 more rows, and 12 more variables: sched_arr_time <int>,
    #> #   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>

    Which is more important: arrival delay or departure delay?

    More context would be helpful to answer this question. Generally departure delay may be considered the most important since it is a proxy for cancelled flights.

  2. Come up with another approach that will give you the same output as not_cancelled %>% count(dest) and not_cancelled %>% count(tailnum, wt = distance) (without using count()).

    #original
    not_cancelled %>% 
      count(dest)
    #> # A tibble: 104 x 2
    #>    dest     n
    #>   <chr> <int>
    #> 1   ABQ   254
    #> 2   ACK   264
    #> 3   ALB   418
    #> 4   ANC     8
    #> 5   ATL 16837
    #> 6   AUS  2411
    #> # ... with 98 more rows
    #another appraoch
    not_cancelled %>%
      group_by(dest) %>%
      summarize(n = n())
    #> # A tibble: 104 x 2
    #>    dest     n
    #>   <chr> <int>
    #> 1   ABQ   254
    #> 2   ACK   264
    #> 3   ALB   418
    #> 4   ANC     8
    #> 5   ATL 16837
    #> 6   AUS  2411
    #> # ... with 98 more rows
    #original
    not_cancelled %>% 
      count(tailnum, wt = distance)
    #> # A tibble: 4,037 x 2
    #>   tailnum      n
    #>     <chr>  <dbl>
    #> 1  D942DN   3418
    #> 2  N0EGMQ 239143
    #> 3  N10156 109664
    #> 4  N102UW  25722
    #> 5  N103US  24619
    #> 6  N104UW  24616
    #> # ... with 4,031 more rows
    #another approach
    not_cancelled %>%
      group_by(tailnum) %>%
      summarize(n = sum(distance))
    #> # A tibble: 4,037 x 2
    #>   tailnum      n
    #>     <chr>  <dbl>
    #> 1  D942DN   3418
    #> 2  N0EGMQ 239143
    #> 3  N10156 109664
    #> 4  N102UW  25722
    #> 5  N103US  24619
    #> 6  N104UW  24616
    #> # ... with 4,031 more rows
  3. Our definition of cancelled flights (is.na(dep_delay) | is.na(arr_delay) ) is slightly suboptimal. Why? Which is the most important column?

    NA’s for dep_delay and dep_time are both the same (8255), whereas arr_time (8713) and arr_delay (9430) likely include missing data that was simply not recorded, or even flights that departed and failed to arrive dur to crashes or other extenuating circumstances. Alternately, for analysis we might want to exlude the arrival na’s in cases where they become contagious.

    summary(flights)
    #>       year          month             day           dep_time   
    #>  Min.   :2013   Min.   : 1.000   Min.   : 1.00   Min.   :   1  
    #>  1st Qu.:2013   1st Qu.: 4.000   1st Qu.: 8.00   1st Qu.: 907  
    #>  Median :2013   Median : 7.000   Median :16.00   Median :1401  
    #>  Mean   :2013   Mean   : 6.549   Mean   :15.71   Mean   :1349  
    #>  3rd Qu.:2013   3rd Qu.:10.000   3rd Qu.:23.00   3rd Qu.:1744  
    #>  Max.   :2013   Max.   :12.000   Max.   :31.00   Max.   :2400  
    #>                                                  NA's   :8255  
    #>  sched_dep_time   dep_delay          arr_time    sched_arr_time
    #>  Min.   : 106   Min.   : -43.00   Min.   :   1   Min.   :   1  
    #>  1st Qu.: 906   1st Qu.:  -5.00   1st Qu.:1104   1st Qu.:1124  
    #>  Median :1359   Median :  -2.00   Median :1535   Median :1556  
    #>  Mean   :1344   Mean   :  12.64   Mean   :1502   Mean   :1536  
    #>  3rd Qu.:1729   3rd Qu.:  11.00   3rd Qu.:1940   3rd Qu.:1945  
    #>  Max.   :2359   Max.   :1301.00   Max.   :2400   Max.   :2359  
    #>                 NA's   :8255      NA's   :8713                 
    #>    arr_delay          carrier              flight       tailnum         
    #>  Min.   : -86.000   Length:336776      Min.   :   1   Length:336776     
    #>  1st Qu.: -17.000   Class :character   1st Qu.: 553   Class :character  
    #>  Median :  -5.000   Mode  :character   Median :1496   Mode  :character  
    #>  Mean   :   6.895                      Mean   :1972                     
    #>  3rd Qu.:  14.000                      3rd Qu.:3465                     
    #>  Max.   :1272.000                      Max.   :8500                     
    #>  NA's   :9430                                                           
    #>     origin              dest              air_time        distance   
    #>  Length:336776      Length:336776      Min.   : 20.0   Min.   :  17  
    #>  Class :character   Class :character   1st Qu.: 82.0   1st Qu.: 502  
    #>  Mode  :character   Mode  :character   Median :129.0   Median : 872  
    #>                                        Mean   :150.7   Mean   :1040  
    #>                                        3rd Qu.:192.0   3rd Qu.:1389  
    #>                                        Max.   :695.0   Max.   :4983  
    #>                                        NA's   :9430                  
    #>       hour           minute        time_hour                  
    #>  Min.   : 1.00   Min.   : 0.00   Min.   :2013-01-01 05:00:00  
    #>  1st Qu.: 9.00   1st Qu.: 8.00   1st Qu.:2013-04-04 13:00:00  
    #>  Median :13.00   Median :29.00   Median :2013-07-03 10:00:00  
    #>  Mean   :13.18   Mean   :26.23   Mean   :2013-07-03 05:02:36  
    #>  3rd Qu.:17.00   3rd Qu.:44.00   3rd Qu.:2013-10-01 07:00:00  
    #>  Max.   :23.00   Max.   :59.00   Max.   :2013-12-31 23:00:00  
    #> 

    Out of curiousity: An alternate way to count na’s by column: https://sebastiansauer.github.io/sum-isna/

    flights %>%
      select(everything()) %>%
      summarize_all(funs(sum(is.na(.))))
    #> # A tibble: 1 x 19
    #>    year month   day dep_time sched_dep_time dep_delay arr_time
    #>   <int> <int> <int>    <int>          <int>     <int>    <int>
    #> 1     0     0     0     8255              0      8255     8713
    #> # ... with 12 more variables: sched_arr_time <int>, arr_delay <int>,
    #> #   carrier <int>, flight <int>, tailnum <int>, origin <int>, dest <int>,
    #> #   air_time <int>, distance <int>, hour <int>, minute <int>,
    #> #   time_hour <int>
  4. Look at the number of cancelled flights per day. Is there a pattern? Is the proportion of cancelled flights related to the average delay?

    Yes, there is a pattern. On a typical day in 2013, the number of flights cancelled was less than .5% and the average dpearture delay was less than 15 minutes.

    However, on both typical days and atypical days, as the percentage of cancellations increases, so does the average delay. This trend holds up to 25% cancellations, whereafter the number of observations becomes exceedingly smaller.

    After 25% cancellations on a given day, one might assume only flights which are not affected extraordinary circumstances are permitted to fly, and therefore are not impacted by excessive delay.

    flights %>%
      group_by(year, month, day) %>%
      summarize(num_flights = n(),
                num_cancelled = sum(is.na(dep_delay)),
                pct_cancelled = mean(is.na(dep_delay)),
                avg_delay     = mean(dep_delay, na.rm = TRUE)
                ) %>%
      ggplot(mapping = aes(x = pct_cancelled * 100, y = avg_delay)) +
      geom_point(aes(alpha = 1/20))+
      geom_smooth(se = FALSE)
    #> `geom_smooth()` using method = 'loess'

    So, out of curiousity, which days had an inordinate number of cancellations, more than 50%?

    flights %>%
      group_by(year, month, day) %>%
      summarize(num_flights = n(),
                num_cancelled = sum(is.na(dep_delay)),
                pct_cancelled = mean(is.na(dep_delay)),
                avg_delay     = mean(dep_delay, na.rm = TRUE)
                ) %>%
      filter(pct_cancelled > .5)
    #> # A tibble: 2 x 7
    #> # Groups:   year, month [1]
    #>    year month   day num_flights num_cancelled pct_cancelled avg_delay
    #>   <int> <int> <int>       <int>         <int>         <dbl>     <dbl>
    #> 1  2013     2     8         930           472     0.5075269  14.85590
    #> 2  2013     2     9         684           393     0.5745614  18.52577

    February 8th and 9th, 2013. What was in the news? Blizzard in the northeast US https://www.theguardian.com/world/2013/feb/08/4000-flights-cancelled-snow-us

  5. Which carrier has the worst delays? Challenge: can you disentangle the effects of bad airports vs. bad carriers? Why/why not? (Hint: think about flights %>% group_by(carrier, dest) %>% summarise(n()))

    not_cancelled %>%
      group_by(carrier) %>%
      summarize(
        avg_dep_delay = mean(dep_delay),
        avg_arr_delay = mean(arr_delay)
      ) %>%
      arrange(desc(avg_arr_delay))
    #> # A tibble: 16 x 3
    #>   carrier avg_dep_delay avg_arr_delay
    #>     <chr>         <dbl>         <dbl>
    #> 1      F9      20.20117      21.92070
    #> 2      FL      18.60598      20.11591
    #> 3      EV      19.83893      15.79643
    #> 4      YV      18.89890      15.55699
    #> 5      OO      12.58621      11.93103
    #> 6      MQ      10.44538      10.77473
    #> # ... with 10 more rows
    
    # Carrier F9 has the worst average arrival delay, being 21.92 minutes
    
    # plot
    plot.data <- not_cancelled %>%
    group_by(carrier) %>%
    summarize(
    avg_dep_delay = mean(dep_delay),
    avg_arr_delay = mean(arr_delay)
    )
    ggplot(plot.data) +
    geom_col(mapping = aes(x = carrier, y = avg_arr_delay))

    
    # effects of bad airports vs. bad carriers
    # we already know carrier F9 has the worst average arrival delay
    
    # which combination of carrier and destination has the worst average arrival delay?
    
    not_cancelled %>%
      group_by(carrier, dest) %>%
      summarize(
        num_flights = n(),
        avg_arr_delay_per_dest = mean(arr_delay)
      ) %>%
      arrange(desc(avg_arr_delay_per_dest))
    #> # A tibble: 312 x 4
    #> # Groups:   carrier [16]
    #>   carrier  dest num_flights avg_arr_delay_per_dest
    #>     <chr> <chr>       <int>                  <dbl>
    #> 1      UA   STL           2              110.00000
    #> 2      OO   ORD           1              107.00000
    #> 3      OO   DTW           2               68.50000
    #> 4      UA   RDU           1               56.00000
    #> 5      EV   CAE         103               42.80583
    #> 6      EV   TYS         313               41.15016
    #> # ... with 306 more rows
    
    # UA into STL, but only 2 flights. 
    
    # let's look at sample of 20 or more flights
    
    not_cancelled %>%
      group_by(carrier, dest) %>%
      summarize(
        num_flights = n(),
        avg_arr_delay_per_dest = mean(arr_delay)
      ) %>%
      arrange(desc(avg_arr_delay_per_dest)) %>%
      filter(num_flights >= 20)
    #> # A tibble: 261 x 4
    #> # Groups:   carrier [16]
    #>   carrier  dest num_flights avg_arr_delay_per_dest
    #>     <chr> <chr>       <int>                  <dbl>
    #> 1      EV   CAE         103               42.80583
    #> 2      EV   TYS         313               41.15016
    #> 3      EV   TUL         294               33.65986
    #> 4      EV   OKC         315               30.61905
    #> 5      MQ   CVG         350               24.02000
    #> 6      EV   MKE        1082               23.23660
    #> # ... with 255 more rows
    
    # EV into CAE tops the list. In fact carrier EV holds the 6 of the top 7 spots.
    # F9 shows up at number 9 into DEN.
    
    # which airport has the worst average arrival delay?
    not_cancelled %>%
      group_by(dest) %>%
      summarize(
        num_flights = n(),
        avg_arr_delay_at_dest = mean(arr_delay)
      ) %>%
      arrange(desc(avg_arr_delay_at_dest))
    #> # A tibble: 104 x 3
    #>    dest num_flights avg_arr_delay_at_dest
    #>   <chr>       <int>                 <dbl>
    #> 1   CAE         106              41.76415
    #> 2   TUL         294              33.65986
    #> 3   OKC         315              30.61905
    #> 4   JAC          21              28.09524
    #> 5   TYS         578              24.06920
    #> 6   MSN         556              20.19604
    #> # ... with 98 more rows
    # Indeed CAE has the worst average arrival delay of any destination
    
    # So what might be causing F9 to top the list of carriers with the worst average arrival delay?
    
    not_cancelled %>%
      group_by(carrier, dest) %>%
      summarize(
        num_flights = n(),
        avg_arr_delay_per_dest = mean(arr_delay)
      ) %>%
      arrange(desc(avg_arr_delay_per_dest)) %>%
      filter(carrier == 'F9')
    #> # A tibble: 1 x 4
    #> # Groups:   carrier [1]
    #>   carrier  dest num_flights avg_arr_delay_per_dest
    #>     <chr> <chr>       <int>                  <dbl>
    #> 1      F9   DEN         681                21.9207
    
    # Carrier F9 only flies into DEN
    
    # What about carrier EV?
    
    not_cancelled %>%
      group_by(carrier, dest) %>%
      summarize(
        num_flights = n(),
        avg_arr_delay_per_dest = mean(arr_delay)
      ) %>%
      arrange(desc(avg_arr_delay_per_dest)) %>%
      filter(carrier == 'EV')
    #> # A tibble: 61 x 4
    #> # Groups:   carrier [1]
    #>   carrier  dest num_flights avg_arr_delay_per_dest
    #>     <chr> <chr>       <int>                  <dbl>
    #> 1      EV   CAE         103               42.80583
    #> 2      EV   TYS         313               41.15016
    #> 3      EV   PBI           6               40.66667
    #> 4      EV   TUL         294               33.65986
    #> 5      EV   OKC         315               30.61905
    #> 6      EV   MKE        1082               23.23660
    #> # ... with 55 more rows
    
    # EV flies into many destinations where their avg arrival delay is relative very high, however it also flies into 51 more destinations. 
    # As yet, we have not compared the delays weighted against number of flights, but let's plot EV's delays per destination.
    
    not_cancelled %>% 
      group_by(carrier, dest) %>%
      summarize(
        num_flights = n(),
        avg_arr_delay_per_dest = mean(arr_delay)
      ) %>%
      filter(carrier == 'EV') %>%
      ggplot(mapping = aes(x = dest, y = avg_arr_delay_per_dest)) +
      geom_col() +
      coord_flip()

    
    # Plot shows most of carrier EV's average arrival delays into its 61 destinations are 20 minutes or less.
    
    # While it has a relatively high average arrival delay as a carrier, and it also holds many top spots 
    
    # For carrier-destination combo's, EV's overall average arrival delay is tempered by the fact it has 61 diverse destinations.
  6. What does the sort argument to count() do. When might you use it?

    Sort = TRUE will sort output in descending order of n

    not_cancelled %>%
      count(dest, sort = TRUE)
    #> # A tibble: 104 x 2
    #>    dest     n
    #>   <chr> <int>
    #> 1   ATL 16837
    #> 2   ORD 16566
    #> 3   LAX 16026
    #> 4   BOS 15022
    #> 5   MCO 13967
    #> 6   CLT 13674
    #> # ... with 98 more rows

    Which is easier than:

    not_cancelled %>%
      group_by(dest) %>%
      summarize(n =n()) %>%
      arrange(desc(n))
    #> # A tibble: 104 x 2
    #>    dest     n
    #>   <chr> <int>
    #> 1   ATL 16837
    #> 2   ORD 16566
    #> 3   LAX 16026
    #> 4   BOS 15022
    #> 5   MCO 13967
    #> 6   CLT 13674
    #> # ... with 98 more rows

    Out of curiosity, what were the maximum arrival delays at each respective destination … ouch:

    not_cancelled %>%
      group_by(dest)%>%
      summarize(
        n = n(),
        max_arr_delay = max(arr_delay)
        )%>%
      arrange(desc(max_arr_delay))
    #> # A tibble: 104 x 3
    #>    dest     n max_arr_delay
    #>   <chr> <int>         <dbl>
    #> 1   HNL   701          1272
    #> 2   CMH  3326          1127
    #> 3   ORD 16566          1109
    #> 4   SFO 13173          1007
    #> 5   CVG  3725           989
    #> 6   TPA  7390           931
    #> # ... with 98 more rows

Grouped mutates (and filters)

Grouping is most useful in conjunction with summarise(), but you can also do convenient operations with mutate() and filter():

  • Find the worst members of each group:

    flights_sml %>% 
      group_by(year, month, day) %>%
      filter(rank(desc(arr_delay)) < 10)
    #> # A tibble: 3,306 x 7
    #> # Groups:   year, month, day [365]
    #>    year month   day dep_delay arr_delay distance air_time
    #>   <int> <int> <int>     <dbl>     <dbl>    <dbl>    <dbl>
    #> 1  2013     1     1       853       851      184       41
    #> 2  2013     1     1       290       338     1134      213
    #> 3  2013     1     1       260       263      266       46
    #> 4  2013     1     1       157       174      213       60
    #> 5  2013     1     1       216       222      708      121
    #> 6  2013     1     1       255       250      589      115
    #> # ... with 3,300 more rows
  • Find all groups bigger than a threshold:

    popular_dests <- flights %>% 
      group_by(dest) %>% 
      filter(n() > 365)
    popular_dests
    #> # A tibble: 332,577 x 19
    #> # Groups:   dest [77]
    #>    year month   day 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      533            529         4      850
    #> 3  2013     1     1      542            540         2      923
    #> 4  2013     1     1      544            545        -1     1004
    #> 5  2013     1     1      554            600        -6      812
    #> 6  2013     1     1      554            558        -4      740
    #> # ... with 3.326e+05 more rows, and 12 more variables:
    #> #   sched_arr_time <int>, 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>
  • Standardize to compute per group metrics:

    popular_dests %>% 
      filter(arr_delay > 0) %>% 
      mutate(prop_delay = arr_delay / sum(arr_delay)) %>% 
      select(year:day, dest, arr_delay, prop_delay)
    #> # A tibble: 131,106 x 6
    #> # Groups:   dest [77]
    #>    year month   day  dest arr_delay   prop_delay
    #>   <int> <int> <int> <chr>     <dbl>        <dbl>
    #> 1  2013     1     1   IAH        11 1.106740e-04
    #> 2  2013     1     1   IAH        20 2.012255e-04
    #> 3  2013     1     1   MIA        33 2.350026e-04
    #> 4  2013     1     1   ORD        12 4.239594e-05
    #> 5  2013     1     1   FLL        19 9.377853e-05
    #> 6  2013     1     1   ORD         8 2.826396e-05
    #> # ... with 1.311e+05 more rows

A grouped filter is a grouped mutate followed by an ungrouped filter. I generally avoid them except for quick and dirty manipulations: otherwise it’s hard to check that you’ve done the manipulation correctly.

Functions that work most naturally in grouped mutates and filters are known as window functions (vs. the summary functions used for summaries). You can learn more about useful window functions in the corresponding vignette: vignette("window-functions").

Exercises

  1. Refer back to the lists of useful mutate and filtering functions. Describe how each operation changes when you combine it with grouping.

    Useful creation functions

    • +, - etc
    • log()
    • lead(), lag()
    • dense_rank(), min_rank(), percent_rank(), row_number(), cume_dist(), ntile()
    • cumsum(), cummean(), cummin(), cummax(), cumany(), cumall()
    • na_if(), coalesce()
    • if_else(), recode(), case_when()
  2. Which plane (tailnum) has the worst on-time record?

    Tailnum N988AT has the worst on-time record in terms of percentage of late arrivals for more than twenty flights.

    not_cancelled %>% 
      group_by(tailnum) %>% 
      summarise(
        n = n(),
        pct_arr_delay = sum(arr_delay > 0) / n()
      ) %>% 
      filter(n > 20) %>% 
      filter(min_rank(desc(pct_arr_delay)) <= 3) %>% 
      arrange(desc(pct_arr_delay))
    #> # A tibble: 3 x 3
    #>   tailnum     n pct_arr_delay
    #>     <chr> <int>         <dbl>
    #> 1  N988AT    35     0.8000000
    #> 2  N983AT    32     0.7500000
    #> 3  N980AT    47     0.7446809

    Alternately, using top_n(1, x), which is the same as if using above filter(min_rank(desc(x)) <= 1)

    not_cancelled %>% 
      group_by(tailnum) %>%
      summarise(
        n= n(),
        late = sum(arr_delay > 0),
        pct_arr_delay =  late / n
        ) %>% 
      filter(n > 20) %>% 
      top_n(1, pct_arr_delay)
    #> # A tibble: 1 x 4
    #>   tailnum     n  late pct_arr_delay
    #>     <chr> <int> <int>         <dbl>
    #> 1  N988AT    35    28           0.8
  3. What time of day should you fly if you want to avoid delays as much as possible?

    On average, avoid 7:00PM to 10:00 PM if you want to avoid delays as much as possible.

    not_cancelled %>% 
      group_by(hour) %>% 
      summarise(
        avg_dep_delay = mean(dep_delay)
      ) %>% 
      filter(min_rank(desc(avg_dep_delay)) <= 4) %>% 
      arrange(desc(avg_dep_delay))
    #> # A tibble: 4 x 2
    #>    hour avg_dep_delay
    #>   <dbl>         <dbl>
    #> 1    19      24.65017
    #> 2    20      24.24992
    #> 3    21      24.16348
    #> 4    17      21.00152
  4. For each destination, compute the total minutes of delay. For each, flight compute the proportion of the total delay for its destination.

    # For each destination, compute the total minutes of delay.
    
    not_cancelled %>% 
      filter(arr_delay > 0) %>% 
      group_by(dest, flight) %>% 
      summarize(
        total_mins = sum(arr_delay)
        #prop_delay = arr_delay / total_mins
      ) %>% 
      select(dest, flight, total_mins) %>% 
      arrange(dest)
    #> # A tibble: 8,505 x 3
    #> # Groups:   dest [103]
    #>    dest flight total_mins
    #>   <chr>  <int>      <dbl>
    #> 1   ABQ     65       1943
    #> 2   ABQ   1505       2544
    #> 3   ACK   1191       1413
    #> 4   ACK   1195         62
    #> 5   ACK   1291        267
    #> 6   ACK   1491       1232
    #> # ... with 8,499 more rows
    by_dest <- not_cancelled %>% 
      filter(arr_delay > 0) %>%   # filter by rows wih arival delays > 0
      group_by(dest) %>%
      mutate(
        total_dest_delay = sum(arr_delay)
        )
    
    # The total_dest_delay gets mutated to each row, respective to the row's dest, and thus is availble to be used below as the detominator in order to calculate the proportional delay.
    
    
    by_flight <- by_dest %>% 
      group_by(flight) %>% 
      mutate(
        prop_delay = sum(arr_delay) / total_dest_delay
        )
    
    # Ungrouping is necessary in order to regroup and summarize. else the results would be a rather uninterptretable as 133K largely non-unique rows, rather than 8K unique rows.  
    
    by_flight %>% 
      ungroup() %>% 
      group_by(dest, flight, total_dest_delay, prop_delay) %>% 
      summarize(
      )
    #> # A tibble: 8,505 x 4
    #> # Groups:   dest, flight, total_dest_delay [?]
    #>    dest flight total_dest_delay prop_delay
    #>   <chr>  <int>            <dbl>      <dbl>
    #> 1   ABQ     65             4487  0.7133943
    #> 2   ABQ   1505             4487  0.6304881
    #> 3   ACK   1191             2974  0.4778077
    #> 4   ACK   1195             2974  0.3190989
    #> 5   ACK   1291             2974  0.1963685
    #> 6   ACK   1491             2974  1.9899126
    #> # ... with 8,499 more rows
  5. Delays are typically temporally correlated: even once the problem that caused the initial delay has been resolved, later flights are delayed to allow earlier flights to leave. Using lag() explore how the delay of a flight is related to the delay of the immediately preceding flight.

    Airports can be big. It may be helpful to have gate information. Lacking gate information, let’s try using as a proxy one carrier at one origin since carriers tend to be assigned one or more specific gates.

    
    not_cancelled %>%
      filter(carrier == "DL", origin == 'LGA', month == 1, dep_delay > 0) %>% 
      arrange(month, day, sched_dep_time) %>% 
      mutate(
        prior_delay = lag(dep_delay)
      ) %>% 
      ggplot() +
      geom_smooth(aes(dep_delay, prior_delay), se = FALSE)
    #> `geom_smooth()` using method = 'loess'

  6. Look at each destination. Can you find flights that are suspiciously fast? (i.e. flights that represent a potential data entry error). Compute the air time a flight relative to the shortest flight to that destination. Which flights were most delayed in the air?

    Using the sigma-three rule as a naive guideline, there are 19 out of 327,000 non-cancelled flights for which the recorded shortest air_time is outside the band of three standard deviations.

    not_cancelled %>%
      group_by(origin, dest) %>% 
      mutate(
        average = mean(air_time, na.rm = TRUE),
        std_dev = sd(air_time),
        x3std_dev = 3 * std_dev,
        avg_less_x3_sd = average - x3std_dev
      ) %>% 
      filter(min_rank(air_time) <= 1, air_time < avg_less_x3_sd) %>% 
      select(origin, dest, average, std_dev, x3std_dev, avg_less_x3_sd, shortest = air_time) %>% 
      arrange(shortest)
    #> # A tibble: 19 x 7
    #> # Groups:   origin, dest [19]
    #>   origin  dest  average  std_dev x3std_dev avg_less_x3_sd shortest
    #>    <chr> <chr>    <dbl>    <dbl>     <dbl>          <dbl>    <dbl>
    #> 1    LGA   BOS 37.85082 4.677608 14.032824       23.81800       21
    #> 2    EWR   SYR 38.98630 2.966448  8.899343       30.08696       30
    #> 3    EWR   BTV 46.25734 3.673901 11.021702       35.23563       34
    #> 4    JFK   ROC 51.88936 4.450083 13.350249       38.53911       35
    #> 5    EWR   RIC 53.54365 5.242607 15.727821       37.81583       35
    #> 6    JFK   ORF 52.58739 4.728371 14.185112       38.40228       36
    #> # ... with 13 more rows
  7. Find all destinations that are flown by at least two carriers. Use that information to rank the carriers.

    # Destinations that are flown by at least two carriers - with summarize()
    
    flights %>%
      group_by(dest) %>% 
      summarize(n_carriers = n_distinct(carrier)) %>% 
      filter(n_carriers >= 2) %>% 
      arrange(dest)
    #> # A tibble: 76 x 2
    #>    dest n_carriers
    #>   <chr>      <int>
    #> 1   ATL          7
    #> 2   AUS          6
    #> 3   AVL          2
    #> 4   BDL          2
    #> 5   BGR          2
    #> 6   BNA          5
    #> # ... with 70 more rows
    # Destinations that are flown by at least two carriers - with mutate(). 
    
    flights %>%
      group_by(dest) %>% 
      mutate(n_carriers = n_distinct(carrier)) %>% 
      filter(rank(n_carriers) >= 2)
    #> # A tibble: 336,774 x 20
    #> # Groups:   dest [103]
    #>    year month   day 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      533            529         4      850
    #> 3  2013     1     1      542            540         2      923
    #> 4  2013     1     1      544            545        -1     1004
    #> 5  2013     1     1      554            600        -6      812
    #> 6  2013     1     1      554            558        -4      740
    #> # ... with 3.368e+05 more rows, and 13 more variables:
    #> #   sched_arr_time <int>, 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>,
    #> #   n_carriers <int>
    # Use that information to rank the carriers. - worst average arrival delay (top 4)
    
    flights %>%
      group_by(dest) %>% 
      mutate(n_carriers = n_distinct(carrier)) %>% 
      filter(n_carriers >= 2) %>% 
      group_by(dest, carrier) %>% 
      summarize(avg_arr_delay = mean(arr_delay,na.rm = TRUE)) %>% 
      select(dest, carrier, avg_arr_delay) %>% 
      filter(min_rank(desc(avg_arr_delay)) <= 4) %>% 
      arrange(dest, desc(avg_arr_delay))
    #> # A tibble: 241 x 3
    #> # Groups:   dest [76]
    #>    dest carrier avg_arr_delay
    #>   <chr>   <chr>         <dbl>
    #> 1   ATL      FL      20.74451
    #> 2   ATL      EV      19.63829
    #> 3   ATL      MQ      14.03400
    #> 4   ATL      UA      10.50000
    #> 5   AUS      AA      16.20950
    #> 6   AUS      B6      11.69272
    #> # ... with 235 more rows
    # Use that information to rank the carriers. - most flights into destination (top 4)
    
    flights %>%
      group_by(dest) %>% 
      mutate(n_carriers = n_distinct(carrier)) %>% 
      filter(n_carriers >= 2) %>% 
      group_by(dest, carrier) %>% 
      summarize(number_of_flights = n()) %>% 
      select(dest, carrier, number_of_flights) %>% 
      filter(min_rank(desc(number_of_flights)) <= 4) %>% 
      arrange(dest, desc(number_of_flights))
    #> # A tibble: 242 x 3
    #> # Groups:   dest [76]
    #>    dest carrier number_of_flights
    #>   <chr>   <chr>             <int>
    #> 1   ATL      DL             10571
    #> 2   ATL      FL              2337
    #> 3   ATL      MQ              2322
    #> 4   ATL      EV              1764
    #> 5   AUS      B6               747
    #> 6   AUS      UA               670
    #> # ... with 236 more rows
  8. For each plane, count the number of flights before the first delay of greater than 1 hour.

    To answer this question, the entire exploratory analysis is left intact to demonstrate the process. Some take-aways include: answer the right question and check your intuitions whenever possible. It would have been easy to stop midway, but something didn’t feel right. Pushing through got the right result.

    The intuition here is to use lead() and probably cumsum() in some way. To do that we will need a temporal sequence of delays based on some other variable. What are the choices for the other variable? Plane equates to tailnum, but does that make sense? Planes arrive at one airport and then on to the next. We’ll see.

    #List the variables.
    
    names(flights)
    #>  [1] "year"           "month"          "day"            "dep_time"      
    #>  [5] "sched_dep_time" "dep_delay"      "arr_time"       "sched_arr_time"
    #>  [9] "arr_delay"      "carrier"        "flight"         "tailnum"       
    #> [13] "origin"         "dest"           "air_time"       "distance"      
    #> [17] "hour"           "minute"         "time_hour"
    # Variables we may want to use:
    # 
    # Year, month and day - for sure
    # sched_dep_time - possibly
    # dep_delay - for sure
    # carrier - likely
    # flight - possibly
    # tailnum - possibly
    # Explore delays by day by plane:
    
    flights %>%
      group_by(year, month, day, tailnum) %>% 
      summarize(total_delay = sum(dep_delay)) %>% 
      filter(month == 7, day == 20) %>%  # filter by an arbitrary day
      arrange(desc(total_delay))
    #> # A tibble: 633 x 5
    #> # Groups:   year, month, day [1]
    #>    year month   day tailnum total_delay
    #>   <int> <int> <int>   <chr>       <dbl>
    #> 1  2013     7    20  N718EV         242
    #> 2  2013     7    20  N121UW         236
    #> 3  2013     7    20  N917XJ         233
    #> 4  2013     7    20  N525MQ         226
    #> 5  2013     7    20  N931XJ         220
    #> 6  2013     7    20  N706JB         210
    #> # ... with 627 more rows
    # On 2018-July-20 tailnum N718EV experienced total departure delays of 242 minutes.
    # Let's explore that further, looking for a pattern we can generalize.
    
    flights %>%
      filter(month == 7, day == 20, tailnum == "N718EV", sched_dep_time, dep_delay) %>% 
      select(tailnum, origin, sched_dep_time, dep_delay, dest, arr_time, arr_delay) %>% 
      arrange(sched_dep_time)
    #> # A tibble: 2 x 7
    #>   tailnum origin sched_dep_time dep_delay  dest arr_time arr_delay
    #>     <chr>  <chr>          <int>     <dbl> <chr>    <int>     <dbl>
    #> 1  N718EV    LGA            930        83   MCI     1251        59
    #> 2  N718EV    LGA           1830       159   RDU     2236       131
    # Maybe tailnum N718EV isn't the best example since it only flew twice that day, both times out of LaGuardia no less. But there may be a pattern here we can use. Let's find a plane with more than two flights in a day and a total departure delay > 60 minutes.
    
    flights %>%
      group_by(year, month, day, tailnum) %>% 
      summarize(
        total_delay = sum(dep_delay),
        n = n()
          ) %>% 
      filter(month == 7, day == 20, n > 2, total_delay > 60) %>%  # arbitrary day
      arrange(desc(total_delay))
    #> # A tibble: 4 x 6
    #> # Groups:   year, month, day [1]
    #>    year month   day tailnum total_delay     n
    #>   <int> <int> <int>   <chr>       <dbl> <int>
    #> 1  2013     7    20  N187JB         121     3
    #> 2  2013     7    20  N591JB         108     3
    #> 3  2013     7    20  N721MQ         102     3
    #> 4  2013     7    20  N374JB          67     4
    # Let's explore tailnum N374JB, which had 4 flights on 2018-July-20 and experienced total departure delays of 67 minutes, looking for a pattern we can generalize.
    
    flights %>%
      filter(month == 7, day == 20, tailnum == "N374JB", sched_dep_time, dep_delay) %>% 
      select(tailnum, origin, sched_dep_time, dep_delay, dest, arr_time, arr_delay) %>% 
      arrange(sched_dep_time)
    #> # A tibble: 4 x 7
    #>   tailnum origin sched_dep_time dep_delay  dest arr_time arr_delay
    #>     <chr>  <chr>          <int>     <dbl> <chr>    <int>     <dbl>
    #> 1  N374JB    JFK            605        -2   MCO      830        -9
    #> 2  N374JB    JFK           1229        10   BUF     1420        29
    #> 3  N374JB    JFK           1635        28   BUF     1817         2
    #> 4  N374JB    JFK           2125        31   BOS     2341        55
    
    # We can see this plane always left out of JFK, but never experienced a departure delay exceeding 60 minutes. So we need to reset on the question: it is not about cumulative delays, but individual delays. 
    # Find a plane which flew multiple times in a day, and experienced at least one departure delay > 60 minutes.
    
    flights %>%
      group_by(year, month, day, tailnum) %>% 
      summarise(
        n_flights = n(),
        max_delay = max(dep_delay)
        ) %>% 
      filter(n_flights > 3, max_delay > 60, month == 7)
    #> # A tibble: 115 x 6
    #> # Groups:   year, month, day [26]
    #>    year month   day tailnum n_flights max_delay
    #>   <int> <int> <int>   <chr>     <int>     <dbl>
    #> 1  2013     7     1  N184JB         4       207
    #> 2  2013     7     1  N273JB         4       124
    #> 3  2013     7     1  N752US         4        81
    #> 4  2013     7     2  N231JB         4       222
    #> 5  2013     7     2  N247JB         4        90
    #> 6  2013     7     2  N281JB         4        87
    #> # ... with 109 more rows
    # On 2013-July-03 tailnum N13988 had 4 flights and experienced a maxiumum departure delay of 140 minutes. Let's explore that one further, again looking for a pattern we can generalize.
    
    flights %>%
      filter(month == 7, day == 03, tailnum == "N13988", sched_dep_time, dep_delay) %>% 
      select(tailnum, origin, sched_dep_time, dep_delay, dest, arr_time, arr_delay) %>% 
      arrange(sched_dep_time)
    #> # A tibble: 4 x 7
    #>   tailnum origin sched_dep_time dep_delay  dest arr_time arr_delay
    #>     <chr>  <chr>          <int>     <dbl> <chr>    <int>     <dbl>
    #> 1  N13988    EWR            650         1   BUF      755       -13
    #> 2  N13988    EWR           1015        92   CHS     1436       147
    #> 3  N13988    EWR           1550       137   CVG     2013       141
    #> 4  N13988    EWR           2055       140   JAX      127       123
    # So, N13988 had one flight before its next dep delay was > 60 mins. Let's try to count that.
    
    flights %>% 
      filter(year == 2013, month == 7, day == 03, tailnum == "N13988") %>% 
      group_by(year, month, day, tailnum) %>% 
      arrange(sched_dep_time) %>%
      mutate(
        flights_before = cumsum(dep_delay <= 60)       # cumsum() - TURNED OUT NOT TO WORK, as below.
        #flights_before = sum(cumall(dep_delay <= 60)) # sum(cumall()) - TURNED out to be KEY to final success
        ) %>% 
      filter(arr_delay > 60 & lag(arr_delay < 60)) %>% 
      select(year, month, day, tailnum, origin, sched_dep_time, dep_delay, flights_before)
    #> # A tibble: 1 x 8
    #> # Groups:   year, month, day, tailnum [1]
    #>    year month   day tailnum origin sched_dep_time dep_delay flights_before
    #>   <int> <int> <int>   <chr>  <chr>          <int>     <dbl>          <int>
    #> 1  2013     7     3  N13988    EWR           1015        92              1
    # That works. Let's try it with one month of data and all planes.
    
    flights %>% 
      filter(year == 2013, month == 7) %>% 
      group_by(year, month, day, tailnum) %>% 
      arrange(sched_dep_time) %>%
      mutate(
        flights_before = cumsum(dep_delay >= 60)
        ) %>% 
      filter(dep_delay > 60 & lag(dep_delay < 60)) %>% 
      select(year, month, day, tailnum, origin, sched_dep_time, dep_delay, flights_before) %>% 
      arrange(desc(flights_before))
    #> # A tibble: 1,086 x 8
    #> # Groups:   year, month, day, tailnum [1,084]
    #>    year month   day tailnum origin sched_dep_time dep_delay flights_before
    #>   <int> <int> <int>   <chr>  <chr>          <int>     <dbl>          <int>
    #> 1  2013     7     9  N738EV    LGA           1500       101              2
    #> 2  2013     7    27  N324JB    JFK           1700       202              2
    #> 3  2013     7     9  N518MQ    JFK           1935        99              2
    #> 4  2013     7     9  N8936A    JFK           2030        99              2
    #> 5  2013     7    11  N179JB    JFK           2030        87              2
    #> 6  2013     7     2  N247JB    JFK           2130        90              2
    #> # ... with 1,080 more rows
    # The results look suspicious in that the greatest number of flights before any plane experienced a departure delay > 60 mins was 2. Let's look at one of those planes.
    
    flights %>%
      filter(month == 7, day == 21, tailnum == "N851MQ", sched_dep_time, dep_delay) %>% 
      select(tailnum, origin, sched_dep_time, dep_delay, dest, arr_time, arr_delay) %>% 
      arrange(sched_dep_time)
    #> # A tibble: 3 x 7
    #>   tailnum origin sched_dep_time dep_delay  dest arr_time arr_delay
    #>     <chr>  <chr>          <int>     <dbl> <chr>    <int>     <dbl>
    #> 1  N851MQ    JFK            815        72   DCA     1025        55
    #> 2  N851MQ    JFK           1720        24   DCA     1936        46
    #> 3  N851MQ    JFK           2140        70   DCA       11        71
    # Having actually checked two planes, they both had 2 flights before experiencing a departure delay > 60 minutes. Intuition is that few if any planes made many more than 3 or 4 flights per day.
    
    flights %>% 
      filter(year == 2013, month == 7) %>% 
      group_by(year, month, day, tailnum) %>% 
      arrange(sched_dep_time) %>%
      mutate(
        flights_before = cumsum(dep_delay >= 60)
        ) %>% 
      filter(dep_delay > 60 & lag(dep_delay < 60)) %>% 
      select(year, month, day, tailnum, origin, sched_dep_time, dep_delay, flights_before) %>% 
      arrange(desc(flights_before))
    #> # A tibble: 1,086 x 8
    #> # Groups:   year, month, day, tailnum [1,084]
    #>    year month   day tailnum origin sched_dep_time dep_delay flights_before
    #>   <int> <int> <int>   <chr>  <chr>          <int>     <dbl>          <int>
    #> 1  2013     7     9  N738EV    LGA           1500       101              2
    #> 2  2013     7    27  N324JB    JFK           1700       202              2
    #> 3  2013     7     9  N518MQ    JFK           1935        99              2
    #> 4  2013     7     9  N8936A    JFK           2030        99              2
    #> 5  2013     7    11  N179JB    JFK           2030        87              2
    #> 6  2013     7     2  N247JB    JFK           2130        90              2
    #> # ... with 1,080 more rows
    # Let's explore the intuition that few planes made more than around 4 flights per day. And if they did, intuition says they are flying short routes, so let's also see what their max_delay was.
    
    flights %>% 
      filter(!is.na(tailnum)) %>% 
      group_by(year, month, day, tailnum) %>% 
      summarize(
        n = n(),
        max_delay = max(dep_delay, na.rm = TRUE)
        ) %>%
      filter(n >= 4) %>% 
      arrange(desc(n))
    #> # A tibble: 2,941 x 6
    #> # Groups:   year, month, day [358]
    #>    year month   day tailnum     n max_delay
    #>   <int> <int> <int>   <chr> <int>     <dbl>
    #> 1  2013     1    16  N15572     5        75
    #> 2  2013     1    28  N26549     5        96
    #> 3  2013     2     2  N329JB     5        78
    #> 4  2013     2     5  N546MQ     5        -1
    #> 5  2013     2     8  N351JB     5        26
    #> 6  2013     2    10  N197JB     5        48
    #> # ... with 2,935 more rows
    # So tailnum N732US made 5 flights on 2013-02-15 and experienced a max departure delay of 114 minutes. Did it make more than 2 flights before this?
    
    flights %>%
      filter(month == 2, day == 15, tailnum == "N732US", sched_dep_time, dep_delay) %>% 
      select(tailnum, origin, sched_dep_time, dep_delay, dest) %>% 
      arrange(sched_dep_time)
    #> # A tibble: 5 x 5
    #>   tailnum origin sched_dep_time dep_delay  dest
    #>     <chr>  <chr>          <int>     <dbl> <chr>
    #> 1  N732US    LGA            700        -4   DCA
    #> 2  N732US    LGA           1100        22   DCA
    #> 3  N732US    LGA           1500        -5   DCA
    #> 4  N732US    LGA           1800        97   DCA
    #> 5  N732US    LGA           2100       114   DCA
    
    # Yes, it made 3 flights before. Therefore, something is wrong with the generalized code above!!!
    # Let's try that one particular plane tailnum N732US made 5 flights on 2013-02-15
    
    flights %>% 
      filter(year == 2013, month == 2, day == 15, tailnum == "N732US") %>% 
      group_by(year, month, day, tailnum) %>% 
      arrange(sched_dep_time) %>%
      mutate(
        flights_before = cumsum(dep_delay <= 60)       # cumsum()
        #flights_before = sum(cumall(dep_delay <= 60)) # sum(cumall()) also works
        ) %>% 
      filter(arr_delay > 60 & lag(arr_delay < 60)) %>% 
      select(year, month, day, tailnum, origin, sched_dep_time, dep_delay, flights_before)
    #> # A tibble: 1 x 8
    #> # Groups:   year, month, day, tailnum [1]
    #>    year month   day tailnum origin sched_dep_time dep_delay flights_before
    #>   <int> <int> <int>   <chr>  <chr>          <int>     <dbl>          <int>
    #> 1  2013     2    15  N732US    LGA           1800        97              3
    
    # That works, which indicates there is in error in the generalized code chunk.
    # We'll try modificating on this generalized chunk, and set the filter for our known observation
    
    flights %>% 
      filter(year == 2013, month == 2, tailnum == "N732US") %>% 
      group_by(year, month, day, tailnum) %>% 
      arrange(sched_dep_time) %>%
      mutate(
      # flights_before = cumsum(lead(dep_delay >= 60)) # DOES NOT WORK
      flights_before = sum(cumall(dep_delay <= 60)) # sum(cumall()) DOES WORK
      # flights_before = cumall(dep_delay <= 60) # DOES NOT WORK
      # flights_before = cumall(lead(dep_delay <= 60)) # DOES NOT WORK
        ) %>% 
      filter(dep_delay > 60 & lag(dep_delay < 60)) %>% 
      select(year, month, day, tailnum, origin, sched_dep_time, dep_delay, flights_before) %>% 
      arrange(desc(flights_before))
    #> # A tibble: 1 x 8
    #> # Groups:   year, month, day, tailnum [1]
    #>    year month   day tailnum origin sched_dep_time dep_delay flights_before
    #>   <int> <int> <int>   <chr>  <chr>          <int>     <dbl>          <int>
    #> 1  2013     2    15  N732US    LGA           1800        97              3
    # Test on larger data set
    
    flights %>% 
      filter(year == 2013, month == 2) %>% 
      group_by(year, month, day, tailnum) %>% 
      arrange(sched_dep_time) %>%
      mutate(
      flights_before = sum(cumall(dep_delay <= 60)) # sum(cumall()) DOES WORK
        ) %>% 
      filter(dep_delay > 60 & lag(dep_delay < 60)) %>% 
      select(year, month, day, tailnum, origin, sched_dep_time, dep_delay, flights_before) %>% 
      arrange(desc(flights_before))
    #> # A tibble: 462 x 8
    #> # Groups:   year, month, day, tailnum [460]
    #>    year month   day tailnum origin sched_dep_time dep_delay flights_before
    #>   <int> <int> <int>   <chr>  <chr>          <int>     <dbl>          <int>
    #> 1  2013     2    28  N12567    EWR           2140        77              4
    #> 2  2013     2    18  N296JB    JFK           1555        94              3
    #> 3  2013     2    15  N732US    LGA           1800        97              3
    #> 4  2013     2    20  N767UW    LGA           1900        93              3
    #> 5  2013     2    21  N756US    LGA           1900        80              3
    #> 6  2013     2    15  N16981    EWR           1940        88              3
    #> # ... with 456 more rows
    
    # Success. Now we're cooking with gas. 
    # Final. All data. Filter results to flights before 60 minute delay to > 3 since this chunk requires relatively more computations than other chunks. 
    
    flights %>% 
      group_by(year, month, day, tailnum) %>% 
      arrange(sched_dep_time) %>%
      mutate(
      flights_before = sum(cumall(dep_delay <= 60)) # sum(cumall()) DOES WORK
        ) %>% 
      filter(dep_delay > 60 & lag(dep_delay < 60) & flights_before > 3) %>% 
      select(year, month, day, tailnum, origin, sched_dep_time, dep_delay, flights_before) %>% 
      arrange(desc(flights_before))
    #> # A tibble: 5 x 8
    #> # Groups:   year, month, day, tailnum [5]
    #>    year month   day tailnum origin sched_dep_time dep_delay flights_before
    #>   <int> <int> <int>   <chr>  <chr>          <int>     <dbl>          <int>
    #> 1  2013     7    29  N329JB    JFK           1901       152              4
    #> 2  2013     7    12  N298JB    JFK           2120        76              4
    #> 3  2013     8    11  N298JB    JFK           2125        72              4
    #> 4  2013     2    28  N12567    EWR           2140        77              4
    #> 5  2013     7     7  N374JB    JFK           2245       227              4