Public Data Dive: 2018 Boeing 737 MAX flights

Public Data Dive: 2018 Boeing 737 MAX flights

The recent grounding of almost all Boeing 737 MAX-series aircraft in the world is, according to a recent CBC commentator, unprecedented. I’m not an aircraft expert (or even a hobbyist), but I do love data and mining publicly-available datasets. Inspired by the nycflights13 R package (a dataset of all the flights in and out of New York City in 2013) and the FlightRadar24 blog post regarding Lion Air flight JT610, I thought I would see what information is accessible to the public about flights that used the 737 MAX-series aircraft.

The nycflights13 R package uses the Bureau of Transportation Statistics Reporting Carrier On-Time Performance and the FAA Releasable Aircraft datasets, both of which are downloadable for free. As per my usual MO, I’ll be using the tidyverse to script up the analysis. Because I’ll be reading in a lot of foreign CSVs with weird names, I’ll also make myself a name cleaner-upper to make nice_columns from BADLY FORMATTED-COLUMNS.

library(tidyverse)
nice_names <- . %>%
  str_to_lower() %>%
  str_replace_all("[- ]", "_")

I opted to get the plane information first. This code was heavily inspired by the scripts in the data-raw/ folder of the nycflights13 R package.

planes_src <- "http://registry.faa.gov/database/yearly/ReleasableAircraft.2018.zip"
planes_tmp <- tempfile(fileext = ".zip")
curl::curl_download(planes_src, planes_tmp)
dir.create("planes")
unzip(planes_tmp, exdir = "planes")

The “planes” folder now contains several oddly named files, of which we are interested in two: “ACFTREF.txt” (airplane model information), and “planes/MASTER.txt” (all registered planes in the US). According to the emergency order by the FAA, the models that were grounded were the Boeing 737-8 and 737-9, although for this analysis I’ll include all 737s for comparison. To do this, I’ll need to build a clean table of plane model info for the planes we’re interested in.

planes_models <- read_csv(
  "planes/ACFTREF.txt",
  col_types = cols(
    .default = col_character(),
    X12 = col_skip()
  ),
  # there are quotation marks as parts of values in this table
  quote = ''
) %>%
  rename_all(nice_names)

planes_models_keep <- planes_models %>%
  filter(mfr == "BOEING", str_detect(model, "^737")) %>%
  mutate(model_label = if_else(str_detect(model, "^737-[89]$"), "737 MAX", "737")) %>%
  select(mfr_mdl_code = code, model_label, model_mfr = mfr, model)

planes_models_keep
## # A tibble: 606 x 4
##    mfr_mdl_code model_label model_mfr model  
##    <chr>        <chr>       <chr>     <chr>  
##  1 1382526      737         BOEING    737-66N
##  2 1382527      737         BOEING    737-7K2
##  3 1382528      737         BOEING    737-9K2
##  4 1382529      737         BOEING    737-8DR
##  5 1382530      737         BOEING    737-8DV
##  6 1382531      737         BOEING    737-7HF
##  7 1382532      737         BOEING    737-46J
##  8 1382540      737         BOEING    737-6Q8
##  9 1383333      737         BOEING    737-2E7
## 10 1384400      737         BOEING    737-7FB
## # … with 596 more rows

Next, I’ll use this to create a table of planes that we’re interested in. I modify the n_number to paste0("N", n_number) here because the tail_num variable in the flights dataset below has the tail number in that form.

planes_master <- read_csv(
  "planes/MASTER.txt",
  col_types = cols(
    .default = col_character(),
    `YEAR MFR` = col_integer(),
    X35 = col_skip()
  )
) %>%
  rename_all(nice_names)

planes_keep <- planes_master %>%
  inner_join(planes_models_keep, by = "mfr_mdl_code") %>%
  mutate(tail_num = paste0("N", n_number)) %>%
  select(tail_num, model_label, model_mfr, model, model_year = year_mfr)

planes_keep
## # A tibble: 2,197 x 5
##    tail_num model_label model_mfr model   model_year
##    <chr>    <chr>       <chr>     <chr>        <int>
##  1 N103WT   737         BOEING    737-700         NA
##  2 N108MS   737         BOEING    737-7BC       2002
##  3 N108WT   737         BOEING    737-700         NA
##  4 N11206   737         BOEING    737-824       2000
##  5 N114JF   737         BOEING    737-4H6       1993
##  6 N12216   737         BOEING    737-824       1998
##  7 N12218   737         BOEING    737-824       1998
##  8 N12221   737         BOEING    737-824       1998
##  9 N12225   737         BOEING    737-824       1998
## 10 N12238   737         BOEING    737-824       1999
## # … with 2,187 more rows

Interestingly, you can use the tail_num to look up flight data from FlightRadar24 like this https://www.flightradar24.com/data/aircraft/<TAIL NUMBER>. You can see an example for one 737 MAX-8 here (note that there aren’t any flights after March 13, when the ban went into effect).

Unfortunately, the pre-zipped option for the on-time performance data was not working at the time of this writing, so I had to download that data by hand (this amounts to one click per month…not too bad for just 2018, but still annoying). I downloaded these files to “flight_zips/2018-{month}.zip”. The code to clean this up was also inspired by the nycflights13 R package cleaning scripts. Note that this is about 300 MB of compressed data and 2.5 GB of uncompressed data before filtering.

extract_month <- function(year = 2018, month) {
  zipfile <- sprintf("flight_zips/%s_%s.zip", year, month)

  files <- unzip(zipfile, list = TRUE)
  # Only extract biggest file
  csv <- files$Name[order(files$Length, decreasing = TRUE)[1]]

  unzip(zipfile, exdir = "flights", junkpaths = TRUE, files = csv)

  src <- paste0("flights/", csv)
  dst <- paste0("flights/", year, "-", month, ".csv")
  file.rename(src, dst)
  dst
}

flight_raw_csvs <- tibble(
  year = 2018,
  month = 1:12,
  csv_file = map2_chr(year, month, extract_month)
)

flight_data <- flight_raw_csvs$csv_file %>%
  map(read_csv, col_types = cols(
    .default = col_character(),
    OP_CARRIER_FL_NUM = col_integer(),
    FL_DATE = col_date(),
    DIVERTED = col_double(),
    AIR_TIME = col_double(),
    CRS_ELAPSED_TIME = col_double(),
    ACTUAL_ELAPSED_TIME = col_double()
  )) %>%
  bind_rows() %>%
  rename_all(nice_names)

flights_keep <- flight_data %>%
  mutate(
    fl_num = paste0(op_unique_carrier, op_carrier_fl_num),
    diverted = if_else(diverted == 1, TRUE, FALSE)
  ) %>%
  select(
    fl_date, fl_num,
    tail_num, origin, dest, diverted,
    schedule_elapsed = crs_elapsed_time,
    actual_elapsed = actual_elapsed_time,
    air_time
  ) %>%
  inner_join(planes_keep, by = "tail_num")

flights_keep

Now that we have a big ol’ data frame of all the flights that took place on a 737 in the US during 2018, we can answer at least one question: were more 737-MAX flights diverted than flights on other 737 aircraft?

flights_keep %>%
  group_by(model_label) %>%
  summarise(
    n = n(),
    n_diverted = sum(diverted),
    pct_diverted = mean(diverted) * 100
  )
## # A tibble: 2 x 4
##   model_label       n n_diverted pct_diverted
##   <chr>         <int>      <int>        <dbl>
## 1 737         2378407       5322        0.224
## 2 737 MAX       38472         98        0.255

The answer appears to be yes, but not by much: 0.25% of 737 MAX flights were diverted, whereas on other 737 flights, 0.22% of flights were diverted. I’m not much of a statistician, but I believe this is a classic case for a Chi-square test, although I think I remember that Chi-squared tests don’t perform well with respect to rare events (flight diversions are definitely rare events). Still, according to this test, even though there were slightly more diverted 737 MAX flights, it is reasonably likely that we would observe this by chance (1 in 5).

observed_diverted <- 98
(observed_not_diverted <- 38472 - 98)
## [1] 38374
(expected_diverted <- 0.2237632 / 100 * 38472)
## [1] 86.08618
(expected_not_diverted <- 38472 - expected_diverted)
## [1] 38385.91
(chi_sq <- sum(
  (observed_diverted - expected_diverted)^2 / expected_diverted,
  (observed_not_diverted - expected_not_diverted)^2 / expected_not_diverted
))
## [1] 1.652501
(p_value <- 1 - pchisq(chi_sq, df = 1))
## [1] 0.1986189

In the off chance this is an artifact of diverted flights being rare, we could also randomly sample the non-MAX flights to see how likely it is that 98 flights of 38472 would be diverted.

non_max_flights <- flights_keep %>%
  filter(model_label == "737") %>%
  select(fl_date, fl_num, diverted)

set.seed(4937)
non_max_flights_boot <- tibble(
  sample_number = 1:1000,
  sample = map(sample_number, ~sample_n(non_max_flights, 38472))
) %>%
  unnest() %>%
  group_by(sample_number) %>%
  summarise(n_diverted = sum(diverted))

ggplot(non_max_flights_boot, aes(n_diverted)) +
  geom_histogram(bins = 30) +
  geom_vline(xintercept = 98) +
  theme_bw()

This is a confirmation of our earlier calculation: it is fairly reasonable that 98 flights of 38472 could be diverted, assuming the two datasets are comparable.

It would be interesting to see similar data as that used by the FlightRadar24 blog post regarding Lion Air flight JT610 (elevation trajectory of the plane during the first 20 minutes of flight) but for all 737 MAX flights. Were there other flights that looked like pilots were fighting with the automatic control system? This data exists at 5-second resolution from FlightRadar24, although even the highest of the paid data plans only lets you download 100 or so of these per day (and there are 2.5 million rows in our flights_keep dataset to examine).

Avatar
Dewey Dunnington
Geoscientist, Programmer, Educator

Related