library(ascii)
library(assertthat)
library(dplyr)
library(ggplot2)
library(knitr)
library(nycflights13)
library(tidyverse)

Note: The stuff in the tidyverse tends to clobber built-in functions, which can make it hard to understand sometimes when things go wrong. If you forget to load dplyer, for instance, there's already a filter function, but it doesn't do the same thing as the one you load with dplyr so you'll get weird error messages if you try to use it. Notably I was getting this message:

Error in filter(airlines, carrier == "US") : object 'carrier' not found

Anyway, now some options for the packages.

options(asciiType="org")
theme_set(theme_minimal(base_family="Palatino"))

And a couple of helpers to print (the tibble objects dump weird characters into the org-mode output).

khead <- function(data) {
  #' Prints a table of the first lines of the dataframe
  #'
  #' @param data The data to print
  table <- kable(head(data))
  table[2] <- gsub(":", "-", table[2])
  table
}

okable <- function(data) {
  #' Prints data as an org-mode table
  #'
  #' @param data The data to print
  table <- kable(data)
  table[2] <- gsub(":", "-", table[2])
  table
}

got_all_your_data <- function(data) {
  #' Checks if there's any missing values
  #'
  #' @param data The vector of data to check
  see_if(noNA(data))
}

filter takes the dataframe as the first argument and every argument after that has to be a boolean expression based on the dataframe. If you pass in more than one condition dplyr will and them. Let's see the flights on Christmas.

christmas <- filter(flights, month==12, day==25)
plot = ggplot(data=christmas) +
  geom_bar(
    mapping=aes(x=origin, fill=carrier)
  ) +
  labs(title="Christmas Flights", x="Origin", y="Count")

ggsave("christmas_flights_counts.png")

Since the default is for the filter to and everything together, if you want to use something else (like OR) you need to create a single argument built with the OR operator (|). If you're using a single column for the OR, there's actually a (possibly) easier way that works like python's in where you check if the value is in a collection. Let's see the flights on Christmas and on Christmas Eve .

both_days <- filter(flights, month==12, day %in% c(24, 25))

plot = ggplot(data=both_days) +
  geom_bar(
    mapping=aes(x=day, fill=origin)
  ) +
  labs(title="Christmas Eve and Day Flights", x="Day", y="Count")

ggsave("christmas_day_and_eve.png")

What if we want New-Years Day and Christmas? In this case we need to keep the month with the day so we can't just use the in operator, but we can use the AND and OR operators to make it work.

and_new_years = filter(flights, (month==12 & day==25) | (month==1 & day==1))

plot = ggplot(data=and_new_years) +
  geom_bar(
    mapping = aes(x=day, fill=origin)
  ) +
  labs(title="Christmas and New Years Day", x="Date", y="Count")

ggsave("christmas_new_years.png")

plot = ggplot(data=and_new_years) +
  geom_boxplot(
    mapping = aes(x=origin, y=dep_delay, color=origin)
  ) +
  labs(title="Christmas and New Years Day Delays", x="Departure Airport", y="Delay (Minutes)")

ggsave("christmas_new_years_delay.png")

What are the shortest flights?

khead(arrange(flights, distance))

The data uses the three-letter airport codes, but there's another dataframe loaded with ggplot2 named airports that we can use to translate them to airport names.

okable(filter(airports, faa %in% c("PHL", "LGA", "EWR")))

So the shortest flights are from Newark Liberty International to La Guardia airport and the Philadelphia International Airport. I wonder why someone would take a plane to go 17 miles? I guess the traffic must be really bad, or some people just have money (and jet fuel) to burn.

To get the data to sort in descending order you have to pass in another function call (desc). Weird.

What are the longest flights?

khead(arrange(flights, desc(distance)))

You can probably guess the airports from the codes, but we'll look them up anyway.

okable(filter(airports, faa %in% c("JFK", "HNL")))

And the longest flights are from John F. Kennedy International airport to Honolulu International airport.

The basic way to get a subset of columns is to pass them in as arguments.

khead(select(flights, month, flight))

If you want to grab a section of contiguous columns you can pass in the first and last ones separated by a colon (the `:` character, not the vital organ).

khead(select(flights, dep_time:arr_time))

Like arrange select uses functions passed in to alter the behavior. Most of them are string-methods-ish.

You can also use logical AND and OR operators or a vector to combine the selections, but we'll just do a basic one-function call. Let's grab all the columns with an underscore in their name.

khead(select(flights, contains("_")))

You can also negate the selections using an exclamation point (!), so we can grab all the columns that don't have an underscore in their name like this.

khead(select(flights, !contains("_")))

That wasn't really helpful, so how about we rename some other columns?

flights <- rename(flights, destination=dest)

This is like the mutate function but it only keeps the new columns you created, not the original columns.

khead(transmute(subset,
       Time_Made_Up=dep_delay - arr_delay,
       Miles_Per_Hour = distance/(air_time/60)
       ))

The reason for summarize is that it works in tandem with groups. Instead of finding the median for the whole column, we can find the median by group.

by_month = group_by(flights, month)

okable(summarize(by_month, median_monthly_delay=median(dep_delay, na.rm=TRUE)))

Is that right? Let's look at a boxplot of the dataset and check.

plot = ggplot(data=by_month) +
  geom_boxplot(
    mapping=aes(x=month, y=dep_delay, fill=factor(month)), 
    outlier.shape=NA,
    na.rm=TRUE
  ) +
  scale_y_continuous(limits = quantile(by_month$dep_delay, c(0.1, 0.9),
                                       na.rm=TRUE)) +
  labs(title="Monthly Delays", x="Month", y="Delay (Minutes)")


ggsave("monthly_delay_boxplot.png")

So, that isn't so exciting, but it does show that the median values tend to be negative (meaning they left early) like our summary said.

dplyr adds an operator called the "pipe" (%>%) that works like a unix pipe does, sending the output of one function to another, allowing you to save some keystrokes adding temporary variables or nesting function calls.

So, here's what we're going to do.

• We're going to pass the flights dataset to the group_by function and group the data by destination
• Then we're going to pass that output to the summarize function and summarize by:
  • the counts
  • mean distance of the flight
  • and the mean arrival delay of each group
• Then we're going to pass the summarized data to filter to get the top 20 that didn't go to Honolulu.

First let's do a check for missing values.

got_all_your_data(flights$arr_delay)
got_all_your_data(flights$distance)

[1] FALSE
attr(,"msg")
[1] "data contains 9430 missing values"
[1] TRUE

There's no missing distances but there are missing arrival delays.

delays_vs_distance <- flights %>%
  group_by(destination) %>%
  summarize(
    count = n(),
    distance = mean(distance, na.rm = TRUE),
    delay = mean(arr_delay, na.rm = TRUE)
  ) %>%
  filter(count > 20, destination != "HNL")

khead(delays_vs_distance)

Note that we're using the mean here, not the median, so there's going to be more variation because of all the outliers. Also, because there are some missing values we had to tell the mean function to ignore them - otherwise the default is to return NA if any values are missing. There weren't any missing distance rows, but passing that argument in doesn't do any harm, I don't think. Now let's plot it.

plot = ggplot(data=delays_vs_distance, mapping=aes(x=distance, y=delay)) +
  geom_point(aes(size=count), alpha=0.4) +
  geom_smooth(se=FALSE) +
  labs(title="Distance Traveled vs Arrival Delay", x="Distance (Miles)",
       y="Delay (Minutes)")

ggsave("distance_vs_delay.png")

It looks like up to a point the further you travel the worse the delay, and then it starts going in the opposite direction, with longer flights having shorter delays in their arrival times.

Note: When you run the first code block to start the session, remember to give it the path to the posts-output folder so you don't need the path when saving the plots.

library("assertthat")
library("ggplot2")
library("knitr")
library("tidyverse")

These are stuff that get re-used in the plotting.

DISPLACEMENT <- "Engine Displacement (Liters)"
EFFICIENCY <-"Gallons Per Mile (Combined)"
DISPLACEMENT_VS_EFFICIENCY = partial(labs, x=DISPLACEMENT, y=EFFICIENCY)

BLUE <- "#4687b7"
RED <- "#ce7b6d"

It's a data frame with 234 rows and 11 variables:

I'll do some column renaming here since I keep forgetting what abbreviations they used.

mpg <- mpg %>%
  rename(
    displacement=displ,
    cylinders=cyl,
    transmission=trans,
    drive_train=drv,
    city_mpg=cty,
    highway_mpg=hwy,
    fuel_type=fl,
    )

str(mpg)

tibble [234 × 11] (S3: tbl_df/tbl/data.frame)
 $ manufacturer: chr [1:234] "audi" "audi" "audi" "audi" ...
 $ model       : chr [1:234] "a4" "a4" "a4" "a4" ...
 $ displacement: num [1:234] 1.8 1.8 2 2 2.8 2.8 3.1 1.8 1.8 2 ...
 $ year        : int [1:234] 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 ...
 $ cylinders   : int [1:234] 4 4 4 4 6 6 6 4 4 4 ...
 $ transmission: chr [1:234] "auto(l5)" "manual(m5)" "manual(m6)" "auto(av)" ...
 $ drive_train : chr [1:234] "f" "f" "f" "f" ...
 $ city_mpg    : int [1:234] 18 21 20 21 16 18 18 18 16 20 ...
 $ highway_mpg : int [1:234] 29 29 31 30 26 26 27 26 25 28 ...
 $ fuel_type   : chr [1:234] "p" "p" "p" "p" ...
 $ class       : chr [1:234] "compact" "compact" "compact" "compact" ...

kable(head(mpg), caption="MPG Dataset (EPA Automobile Data)")

Table: MPG Dataset (EPA Automobile Data)

The EPA uses a weighted average to create a combined metric that uses both Highway and City mileage.

\[ \textrm{Combined Fuel Economy} = 0.55 \times \textit{city} + 0.45 \times \textit{highway} \]

I've also seen it done as the harmonic mean (on Illustrative Mathematics), but if you round to the nearest whole number it comes out about the same.

mpg$combined_fuel_economy <- 0.55 * mpg$city_mpg + 0.45 * mpg$highway_mpg

plot = ggplot(data=mpg) +
  geom_point(mapping=aes(x=displacement, y=combined_fuel_economy, color=combined_fuel_economy)) +
  labs(x=DISPLACEMENT, y="Combined MPG",
       title="Displacement vs Mileage")

filename = "displacement_vs_combined_mileage.png"
ggsave(filename)

We can plot them together to see how much the combined MPG differs from the highway MPG.

plot = ggplot(data=mpg) +
  geom_smooth(mapping=aes(x=displacement, y=combined_fuel_economy), color="red") +
  geom_smooth(mapping=aes(x=displacement, y=highway_mpg), color="blue") +
  labs(x=DISPLACEMENT, y="MPG",
       title="Displacement vs Mileage (Highway and Combined)")

filename = "highway_vs_combined_mileage.png"
ggsave(filename)

Adding the city mileage drops the mileage somewhat, although the line looks to be about the same shape.

The EPA's explanation of their vehicle labels mentions that the reason why they put a Fuel Consumption Rate is that the Miles Per Gallon metric is actually not a good way to compare mileage - it is subject to a "Miles Per Gallon Illusion" (ReaClearScience has an article about it). Here's how the number of gallons used per 1,000 miles driven changes with miles-per-gallon.

mileage <- seq(1, 100)
miles = 1000
gallons_per_1000_miles <- miles/mileage

improvement <- data.frame(MPG=mileage,
                          gallons=gallons_per_1000_miles)
plot = ggplot(data=improvement) +
  geom_point(mapping=aes(x=mileage, y=gallons), alpha=0.5, color="blue") +
  geom_line(mapping=aes(x=mileage, y=gallons), color="blue") +
  labs(title="Fuel Consumption vs MPG",
       x="Miles Per Gallon",
       y="Gallons Used Per 1,000 Miles")

ggsave("gallons_vs_mpg.png")

As you can see the improvement in the number of gallons used goes down pretty fast as the miles-per-gallon goes up. Let's see about adding a fuel-efficency column. What happens when you switch to Gallons Per Mile?

gallons_per_mile <- mileage
gallons_per_1000_miles <- miles * gallons_per_mile

improvement <- data.frame(gpm=gallons_per_mile,
                          gallons=gallons_per_1000_miles)
plot = ggplot(data=improvement) +
  geom_point(mapping=aes(x=gpm, y=gallons), alpha=0.5, color="blue") +
  geom_line(mapping=aes(x=gpm, y=gallons), color="blue") +
  labs(title="Fuel Consumption vs Gallons Per Mile",
       x="Gallons Per Mile",
       y="Gallons Used Per 1,000 Miles")

ggsave("gallons_vs_efficiency.png")

I suppose you could just figure that that'd be the case, but this is about plotting.

mpg$gallons_per_mile <- 1/mpg$combined_fuel_economy

plot = ggplot(data=mpg, mapping=aes(x=displacement, y=gallons_per_mile)) +
  geom_point() +
  DISPLACEMENT_VS_EFFICIENCY(title="Displacement vs Fuel Efficiency")

filename = "displacement_vs_fuel_efficiency.png"
ggsave(filename)

Even though this is just the inverse of the MPG plot it seems more intuitive to me. I didn't really notice that big outlier at the top when I plotted it using MPG, so I think I'll stick with this.

plot = ggplot(data=mpg) +
  geom_point(mapping=aes(x=displacement, y=gallons_per_mile, color=class),
             position="jitter") +
  DISPLACEMENT_VS_EFFICIENCY(title="Displacement vs Fuel Efficiency")

filename = "displacement_vs_mileage_with_class.png"
ggsave(filename)

I added a little jitter because it looked like the top outlier was just one point. It looks like 2-seaters and one of the mid-sized cars had unusually good mileage given their engine displacement. What's a two-seater?

kable(filter(mpg, class=="2seater"))

Just Corvettes here.

While the use of color is useful for comparing the classes on the same plot, sometimes it can be easier to interpret plots that have the categorical classes separated into sub-plots (called facets).

By default the faceting functions use aphabetical order, I'm going to use a slightly different ordering to see if it makes it clearer how the classes compare (this is based on the medians that I plotted below).

mpg$class_order <- factor(mpg$class,
                          levels=c("compact", "subcompact", "midsize",
                                   "2seater",
                                   "minivan", "suv", "pickup"))

plot = ggplot(data=mpg) +
  geom_point(mapping=aes(x=displacement, y=gallons_per_mile, color=gallons_per_mile)) +
  facet_grid(. ~ class_order) + 
  DISPLACEMENT_VS_EFFICIENCY(title="Displacement vs Efficiency")

filename = "displacement_vs_mileage_faceted_classes.png"
ggsave(filename)

I'm using facet_grid here, which is actually meant to use two variables (creating row and column separation). By putting the . before the ~ we're telling it to only facet using columns. There's also a facet_wrap function that only does one variable but adds the ability to specify the number of rows you want to use, which might be helpful if you have a lot of classes.

I wouldn't have thought it, but the Corvette falls between the minivan and the sub-compact for fuel efficiency. I guess when you consider that it can only be used to transport two people and not anything relatively large it's still inefficient, just not as much if you're only using it to commute.

Instead of plotting a scatterpoint, you can use geom_smooth to fit a single line representing the trend for the data (using a linear model).

plot = ggplot(data=mpg, mapping=aes(x=displacement, y=gallons_per_mile)) +
  geom_smooth() +
  geom_point() +
  facet_grid(class_order ~ drive_train) +
  DISPLACEMENT_VS_EFFICIENCY(title="Displacement vs Efficiency By Type")

ggsave("displacement_vs_mileage_smooth.png")

One interesting thing here is that the four-wheel-drive compact and sub-compact cars did better than some of the front-wheel drive cars in the same class.

If you look at the front-wheel-drive sub-compact cars you'll see that the curve-fitting can get a little funky if you're not careful.

Note the syntax for the plot. Before we were passing the aes mapping to the geometry object. You can do that here too, but if you pass it to the ggplot call then it will automatically apply it to all the geometries. If you wanted to then change one of them you could pass in the aes mapping to it and override the base mapping you gave to ggplot.

Besides passing in different aes mappings, you can also pass in different sub-sets of the data to each geometry (using dplyr's filter to isolate the SUVs here).

plot = ggplot(data=mpg, mapping=aes(x=displacement, y=gallons_per_mile)) +
  geom_smooth(data=filter(mpg, class=="suv")) +
  geom_point(mapping=aes(color=class), position="jitter") +
  DISPLACEMENT_VS_EFFICIENCY(title="Displacement vs Fuel Efficiency")

filename = "displacement_vs_mileage_subsets.png"
ggsave(filename)

Looking at the trend-line for SUVs you can see that it encompasses a surprising swath of the displacements, and it is mostly upwardly linear to begin with but then flattens out.

plot = ggplot(
  data = mpg,
  mapping = aes(x = displacement, y = gallons_per_mile, color = drive_train)
) +
  geom_point() +
  geom_smooth(se = FALSE) +
  DISPLACEMENT_VS_EFFICIENCY(title="Effect of Drive (Four, Front, or Rear) on Mileage")

ggsave("displacement_exercise_2.png")

It kind of looks like an emerging pattern that up to a point gasoline consumption goes up with engine size but then once a certain point is passed consumption flattens out, even as the engines get bigger.

Just out of curiosity, I'd like to see what those small four-wheel drive cars are.

small <- filter(mpg, displacement<3 & drive_train=="Four-Wheel")
plot = ggplot(
  data=small,
  mapping=aes(x=displacement, y=gallons_per_mile, color=model)
) +
  geom_point() +
  facet_grid(class_order ~ manufacturer) +
  DISPLACEMENT_VS_EFFICIENCY(title="4WD Less Than 3 Liters")

ggsave("small_4wd.png")

There are only three companies making these small four-wheel-drive automobiles, with Audi being the one focused on sedans.

The Bar Plot is an example of what the Grammar of Graphics calls a Statistical Transformation (I guess the scatter plots are so called "Identity" transformations, but, whatever). It doesn't display the values in the data (unless you tell it to) but instead shows how often categorical values appear in the data. We can also have it set the color to create a stacked bar plot.

plot = ggplot(data=mpg) +
  geom_bar(
    mapping = aes(x=class_order, fill=drive_train)
  ) +
  labs(title="Count of Cars by Class", x="Class Ordered by Median MPG", y="Count")

ggsave("class_count.png")

So most models of car are SUVs, how disheartening. This is only a sub-set of cars, though, and it doesn't reflect sales, just model counts.

Instead of using counts we can use proportions to show how the sub-category portions compare across categories.

plot = ggplot(data=mpg) +
  geom_bar(
    mapping = aes(x=class_order, fill=drive_train),
    position="fill"
  ) +
  labs(title="Proportion of Drives by Class", x="Class Ordered by Median MPG", y="Proportion")

ggsave("class_count_filled.png")

This shows what proportion if each class of car has what type of drive.

To make it easier to see how sub-categories compare within the class, we can change the positioning so that the bars are put side-by-side.

plot = ggplot(data=mpg) +
  geom_bar(
    mapping = aes(x=class_order, fill=drive_train),
    position="dodge"
  ) +
  labs(title="Count of Model Drives by Class", x="Class Ordered by Median MPG", y="Count")

ggsave("class_count_dodged.png")

Our original bar graph showed the counts, but we can show the same plot using the proportion that each categories' count is of the total count as the y-axis.

plot = ggplot(data=mpg, mapping=aes(color=class)) +
  geom_bar(
    mapping = aes(x=class_order, y=..prop.., group=1)
  ) +
  labs(title="Proportion of Cars by Class", x="Class Ordered by Median MPG", y="Proportion")

ggsave("class_proportion.png")

Two things:

• That funky ..prop.. is really how R specifies that argument
• The group=1 argument tells ggplot that there's only one group for the whole dataset, otherwise each category is taken as a separate group and they are all end up as 100% of their group.
• I couldn't get ggplot to set the colors based on the class

That was a nice clean way to show it, and it's flexible, but maybe a good old-fashioned box-plot is the best, in the end.

plot = ggplot(data=mpg) +
  geom_boxplot(
    mapping = aes(x=class_order, y=gallons_per_mile, color=class_order),
  ) +
  labs(title="Spread of Efficiency by Class", y=EFFICIENCY, x="Class")

ggsave("mpg_class_boxplot.png")

library(assertthat)
library(knitr)
library(lubridate)
library(tidyverse)

Lubridate is part of the tidyverse but it has to be loaded for some reason.

DATA_URL <- "https://data.cityofnewyork.us/api/views/833y-fsy8/rows.csv?accessType=DOWNLOAD" 

data = read_csv(DATA_URL)

[1mindexed[0m [32m0B[0m in [36m 0s[0m, [32m0B/s[0m
[1mindexed[0m [32m393.21kB[0m in [36m 0s[0m, [32m983.15kB/s[0m
                                                                              
[1mindexed[0m [32m524.28kB[0m in [36m 1s[0m, [32m1.03MB/s[0m
[1mindexed[0m [32m655.36kB[0m in [36m 1s[0m, [32m1.22MB/s[0m
[1mindexed[0m [32m786.43kB[0m in [36m 1s[0m, [32m1.24MB/s[0m
[1mindexed[0m [32m917.50kB[0m in [36m 1s[0m, [32m1.23MB/s[0m
                                                                              
[1mindexed[0m [32m1.05MB[0m in [36m 1s[0m, [32m1.18MB/s[0m
[1mindexed[0m [32m1.18MB[0m in [36m 1s[0m, [32m1.18MB/s[0m
[1mindexed[0m [32m1.31MB[0m in [36m 1s[0m, [32m1.06MB/s[0m
[1mindexed[0m [32m1.44MB[0m in [36m 2s[0m, [32m947.54kB/s[0m
                                                                              
[1mindexed[0m [32m1.57MB[0m in [36m 2s[0m, [32m1.00MB/s[0m
[1mindexed[0m [32m1.70MB[0m in [36m 2s[0m, [32m1.04MB/s[0m
[1mindexed[0m [32m1.83MB[0m in [36m 2s[0m, [32m1.04MB/s[0m
[1mindexed[0m [32m1.97MB[0m in [36m 2s[0m, [32m1.02MB/s[0m
[1mindexed[0m [32m2.10MB[0m in [36m 2s[0m, [32m1.03MB/s[0m
[1mindexed[0m [32m2.23MB[0m in [36m 2s[0m, [32m994.69kB/s[0m
[1mindexed[0m [32m2.36MB[0m in [36m 2s[0m, [32m965.18kB/s[0m
[1mindexed[0m [32m2.49MB[0m in [36m 3s[0m, [32m939.81kB/s[0m
[1mindexed[0m [32m2.62MB[0m in [36m 3s[0m, [32m951.26kB/s[0m
[1mindexed[0m [32m2.75MB[0m in [36m 3s[0m, [32m919.65kB/s[0m
[1mindexed[0m [32m2.88MB[0m in [36m 3s[0m, [32m855.66kB/s[0m
[1mindexed[0m [32m3.01MB[0m in [36m 3s[0m, [32m892.99kB/s[0m
[1mindexed[0m [32m3.15MB[0m in [36m 3s[0m, [32m904.11kB/s[0m
[1mindexed[0m [32m3.28MB[0m in [36m 4s[0m, [32m867.67kB/s[0m
[1mindexed[0m [32m3.41MB[0m in [36m 4s[0m, [32m878.71kB/s[0m
[1mindexed[0m [32m3.54MB[0m in [36m 4s[0m, [32m888.47kB/s[0m
[1mindexed[0m [32m3.67MB[0m in [36m 4s[0m, [32m899.43kB/s[0m
[1mindexed[0m [32m3.80MB[0m in [36m 4s[0m, [32m846.12kB/s[0m
[1mindexed[0m [32m3.93MB[0m in [36m 5s[0m, [32m861.49kB/s[0m
[1mindexed[0m [32m4.06MB[0m in [36m 5s[0m, [32m876.08kB/s[0m
[1mindexed[0m [32m4.19MB[0m in [36m 5s[0m, [32m854.40kB/s[0m
[1mindexed[0m [32m4.33MB[0m in [36m 5s[0m, [32m877.77kB/s[0m
[1mindexed[0m [32m4.46MB[0m in [36m 5s[0m, [32m855.56kB/s[0m
[1mindexed[0m [32m4.59MB[0m in [36m 5s[0m, [32m865.02kB/s[0m
[1mindexed[0m [32m4.72MB[0m in [36m 6s[0m, [32m855.32kB/s[0m
[1mindexed[0m [32m4.85MB[0m in [36m 6s[0m, [32m847.74kB/s[0m
[1mindexed[0m [32m4.98MB[0m in [36m 6s[0m, [32m870.08kB/s[0m
[1mindexed[0m [32m5.11MB[0m in [36m 6s[0m, [32m892.61kB/s[0m
[1mindexed[0m [32m5.11MB[0m in [36m 6s[0m, [32m892.93kB/s[0m
[1mindexed[0m [32m1.00TB[0m in [36m 6s[0m, [32m174.56GB/s[0m
                                                                              
[1mRows: [22m[34m23585[39m [1mColumns: [22m[34m19[39m
[36m──[39m [1mColumn specification[22m [36m──────────────────────────────────────────────────────────────────────────────[39m
[1mDelimiter:[22m ","
[31mchr[39m  (10): OCCUR_DATE, BORO, LOCATION_DESC, PERP_AGE_GROUP, PERP_SEX, PERP_R...
[32mdbl[39m   (7): INCIDENT_KEY, PRECINCT, JURISDICTION_CODE, X_COORD_CD, Y_COORD_CD...
[33mlgl[39m   (1): STATISTICAL_MURDER_FLAG
[34mtime[39m  (1): OCCUR_TIME

[36mℹ[39m Use [30m[47m[30m[47m`spec()`[47m[30m[49m[39m to retrieve the full column specification for this data.
[36mℹ[39m Specify the column types or set [30m[47m[30m[47m`show_col_types = FALSE`[47m[30m[49m[39m to quiet this message.

spec(data)

cols(
  INCIDENT_KEY = [32mcol_double()[39m,
  OCCUR_DATE = [31mcol_character()[39m,
  OCCUR_TIME = [34mcol_time(format = "")[39m,
  BORO = [31mcol_character()[39m,
  PRECINCT = [32mcol_double()[39m,
  JURISDICTION_CODE = [32mcol_double()[39m,
  LOCATION_DESC = [31mcol_character()[39m,
  STATISTICAL_MURDER_FLAG = [33mcol_logical()[39m,
  PERP_AGE_GROUP = [31mcol_character()[39m,
  PERP_SEX = [31mcol_character()[39m,
  PERP_RACE = [31mcol_character()[39m,
  VIC_AGE_GROUP = [31mcol_character()[39m,
  VIC_SEX = [31mcol_character()[39m,
  VIC_RACE = [31mcol_character()[39m,
  X_COORD_CD = [32mcol_double()[39m,
  Y_COORD_CD = [32mcol_double()[39m,
  Latitude = [32mcol_double()[39m,
  Longitude = [32mcol_double()[39m,
  Lon_Lat = [31mcol_character()[39m
)

floor_date is a function from lubridate that comes as part of the tidyverse. It's called floor because it acts like the floor function in math (round down).

Note: For my future self: using the dataclass is the python way to make a data object, but numba doesn't like it, so we can't use it.

@dataclass
class InsertionOutput:
    element_count: int
    comparisons: int
    swaps: int
    elements: MutableSequence

I'll use the regular namedtuple instead.

InsertionOutput = namedtuple("InsertionOutput", ["element_count",
                                                 "comparisons",
                                                 "swaps",
                                                 "elements"])

def insertion_sort(elements: MutableSequence) -> InsertionOutput:
    """Sorts elements using iterative insertion-sort

    Args:
     elements: sortable collection of elements

    Returns:
     count of elements, comparisons made, swaps made, sorted elements
    """
    comparisons = swaps = 0
    for next_unsorted_cell in range(1, len(elements)):
        thing_to_insert = elements[next_unsorted_cell]

        in_front_of_me, to_the_right = (next_unsorted_cell - 1,
                                        next_unsorted_cell)

        while not (in_front_of_me < 0 or
                   elements[in_front_of_me] <= thing_to_insert):
            comparisons += 1
            swaps += 1
            elements[to_the_right] = elements[in_front_of_me]
            in_front_of_me, to_the_right = (in_front_of_me - 1,
                                            in_front_of_me)

        elements[to_the_right] = thing_to_insert
        swaps += 1

    return InsertionOutput(len(elements), comparisons, swaps, elements)

I negated the while-condition and re-stated the body to make more sense to me. Hopefully it's still clear what's going on.

numba_sort = njit(insertion_sort)
things_to_sort = [numba_random.integers(low=0, high=count, size=count)
                  for count in range(1, 10**5 + 1, 1000)]

with TIMER:
    sort_output = Parallel(n_jobs=-1)(
        delayed(numba_sort)(thing_to_sort)
        for thing_to_sort in things_to_sort)

Started: 2022-01-16 22:55:12.336777
Ended: 2022-01-16 22:55:23.534382
Elapsed: 0:00:11.197605

SIZE, COMPARISONS, SWAPS = 0, 1, 2
unzipped = list(zip(*sort_output))
count_frame = pandas.DataFrame({"Elements": unzipped[SIZE],
                                "Insertion": unzipped[SWAPS]})

count_frame["n^2"] = count_frame.Elements**2

melted = count_frame.melt(id_vars=["Elements"], var_name="Source", value_name="Runtime")