rud.is

Over on The DO Loop, @RickWicklin does a nice job [visualizing the causes of airline crashes](http://blogs.sas.com/content/iml/2015/03/30/visualizing-airline-crashes/) in SAS using a mosaic plot. More often than not, I find mosaic plots can be a bit difficult to grok, but Rick’s use was spot on and I believe it shows the data pretty well, but I also thought I’d take the opportunity to:

– Give @jennybc’s new [googlesheets](http://github.com/jennybc/googlesheets) a spin
– Show some `dplyr` & `tidyr` data wrangling (never can have too many examples)
– Crank out some `ggplot` zero-based streamgraph-y area charts for the data with some extra `ggplot` wrangling for good measure

I also decided to use the colors in the [original David McCandless/Kashan visualization](http://www.informationisbeautiful.net/visualizations/plane-truth-every-single-commercial-plane-crash-visualized/).

#### Getting The Data

As I mentioned, @jennybc made a really nice package to interface with Google Sheets, and the IIB site [makes the data available](https://docs.google.com/spreadsheet/ccc?key=0AjOUPqcIwvnjdEx2akx5ZjJXSk9oM1E3dWpqZFJ6Nmc&usp=drive_web#gid=1), so I copied it to my Google Drive and gave her package a go:

library(googlesheets)
library(ggplot2) # we'll need the rest of the libraries later
library(dplyr)   # but just getting them out of the way
library(tidyr)
 
# this will prompt for authentication the first time
my_sheets <- list_sheets()
 
# which one is the flight data one
grep("Flight", my_sheets$sheet_title, value=TRUE)
 
## [1] "Copy of Flight Risk JSON" "Flight Risk JSON" 
 
# get the sheet reference then the data from the second tab
flights <- register_ss("Flight Risk JSON")
flights_csv <- flights %>% get_via_csv(ws = "93-2014 FINAL")
 
# take a quick look
glimpse(flights_csv)
 
## Observations: 440
## Variables:
## $ date       (chr) "d", "1993-01-06", "1993-01-09", "1993-01-31", "1993-02-08", "1993-02-28", "...
## $ plane_type (chr) "t", "Dash 8-311", "Hawker Siddeley HS-748-234 Srs", "Shorts SC.7 Skyvan 3-1...
## $ loc        (chr) "l", "near Paris Charles de Gualle", "near Surabaya Airport", "Mt. Kapur", "...
## $ country    (chr) "c", "France", "Indonesia", "Indonesia", "Iran", "Taiwan", "Macedonia", "Nor...
## $ ref        (chr) "r", "D-BEAT", "PK-IHE", "9M-PID", "EP-ITD", "B-12238", "PH-KXL", "LN-TSA", ...
## $ airline    (chr) "o", "Lufthansa Cityline", "Bouraq Indonesia", "Pan Malaysian Air Transport"...
## $ fat        (chr) "f", "4", "15", "14", "131", "6", "83", "3", "6", "2", "32", "55", "132", "4...
## $ px         (chr) "px", "20", "29", "29", "67", "22", "56", "19", "22", "17", "38", "47", "67"...
## $ cat        (chr) "cat", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A2", "A1", "A1", "A1...
## $ phase      (chr) "p", "approach", "initial_climb", "en_route", "en_route", "approach", "initi...
## $ cert       (chr) "cert", "confirmed", "probable", "probable", "confirmed", "probable", "confi...
## $ meta       (chr) "meta", "human_error", "mechanical", "weather", "human_error", "weather", "h...
## $ cause      (chr) "cause", "pilot & ATC error", "engine failure", "low visibility", "pilot err...
## $ notes      (chr) "n", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
 
# the spreadsheet has a "helper" row for javascript, so we nix it
flights_csv <- flights_csv[-1,] # js vars removal
 
# and we convert some columns while we're at it
flights_csv %>%
  mutate(date=as.Date(date),
         fat=as.numeric(fat),
         px=as.numeric(px)) -> flights_csv

#### A Bit of Cleanup

Despite being a spreadsheet, the data needs some cleanup and there’s no real need to include “grounded” or “unknown” in the flight phase given the limited number of incidents in those categories. I’d actually mention that descriptively near the visual if this were anything but a blog post.

The area chart also needs full values for each category combo per year, so we use `expand` from `tidyr` with `left_join` and `mutate` to fill in the gaps.

Finally, we make proper, ordered labels:

flights_csv %>%
  mutate(year=as.numeric(format(date, "%Y"))) %>%
  mutate(phase=tolower(phase),
         phase=ifelse(grepl("take", phase), "takeoff", phase),
         phase=ifelse(grepl("climb", phase), "takeoff", phase),
         phase=ifelse(grepl("ap", phase), "approach", phase)) %>%
  count(year, meta, phase) %>%
  left_join(expand(., year, meta, phase), ., c("year", "meta", "phase")) %>% 
  mutate(n=ifelse(is.na(n), 0, n)) %>% 
  filter(!phase %in% c("grounded", "unknown")) %>%
  mutate(phase=factor(phase, 
                      levels=c("takeoff", "en_route", "approach", "landing"),
                      labels=c("Takeoff", "En Route", "Approach", "Landing"),
                      ordered=TRUE)) -> flights_dat

I probably took some liberties lumping “climb” in with “takeoff”, but I’d’ve asked an expert for a production piece just as I would hope folks doing work for infosec reports or visualizations would consult someone knowledgable in cybersecurity.

#### The Final Plot

I’m a big fan of an incremental, additive build idiom for `ggplot` graphics. By using the `gg <- gg + …` style one can move lines around, comment them out, etc without dealing with errant `+` signs. It also forces a logical separation of ggplot elements. Personally, I tend to keep my build orders as follows: - main `ggplot` call with mappings if the graph is short, otherwise add the mappings to the `geom`s - all `geom_` or `stat_` layers in the order I want them, and using line breaks to logically separate elements (like `aes`) or to wrap long lines for easier readability. - all `scale_` elements in order from axes to line to shape to color to fill to alpha; I'm not as consistent as I'd like here, but keeping to this makes it really easy to quickly hone in on areas that need tweaking - `facet` call (if any) - label setting, always with `labs` unless I really have a need for using `ggtitle` - base `theme_` call - all other `theme` elements, one per `gg <- gg +` line I know that's not everyone's cup of tea, but it's just how I roll `ggplot`-style. For this plot, I use a smoothed stacked plot with a custom smoother and also use Futura Medium for the text font. Substitute your own fav font if you don't have Futura Medium.

flights_palette <- c("#702023", "#A34296", "#B06F31", "#939598", "#3297B0")
 
gg <- ggplot(flights_dat, aes(x=year, y=n, group=meta)) 
gg <- gg + stat_smooth(mapping=aes(fill=meta), geom="area",
                       position="stack", method="gam", formula=y~s(x)) 
gg <- gg + scale_fill_manual(name="Reason:", values=flights_palette, 
                             labels=c("Criminal", "Human Error",
                                      "Mechanical", "Unknown", "Weather"))
gg <- gg + scale_y_continuous(breaks=c(0, 5, 10, 13))
gg <- gg + facet_grid(~phase)
gg <- gg + labs(x=NULL, y=NULL, title="Crashes by year, by reason & flight phase")
gg <- gg + theme_bw()
gg <- gg + theme(legend.position="bottom")
gg <- gg + theme(text=element_text(family="Futura Medium"))
gg <- gg + theme(plot.title=element_text(face="bold", hjust=0))
gg <- gg + theme(panel.grid=element_blank())
gg <- gg + theme(panel.border=element_blank())
gg <- gg + theme(strip.background=element_rect(fill="#525252"))
gg <- gg + theme(strip.text=element_text(color="white"))
gg

That ultimately produces:

with the facets ordered by takeoff, flying, approaching landing and actual landing phases. Overall, things have gotten way better, though I haven’t had time to look in to the _bump_ between 2005 and 2010 for landing crashes.

As an aside, Boeing has a [really nice PDF](http://www.boeing.com/news/techissues/pdf/statsum.pdf) on some of this data with quite a bit more detail.

Vis expert Naomi Robbins did an excellent [critique](http://www.forbes.com/sites/naomirobbins/2015/03/19/color-problems-with-figures-from-the-jerusalem-post/) of the [graphics](http://www.jpost.com/Israel-Elections/Analysis-The-Israel-election-decided-by-one-vote-394229) that went along with an article on Israeli election in the Jerusalem Post.

Non-uniform and color-blind-unfriendly categorical colors and disproportionate arc sizes are definitely three substantial issues in that series of visualizations. We can rectify all of them with two new packages of mine: [waffle](http://github.com/hrbrmstr/waffle) & [adobecolor](http://github.com/hrbrmstr/adobecolor). The former provides a good alternative to pie charts (no charts at all are a good alternative to pie charts) and the latter makes it possible to share color palettes without passing long strings of hex-encoded colors.

Using [XScope](http://xscopeapp.com/) I encoded a color-blind-friendly palette from [Brian Connelly](http://bconnelly.net/2013/10/creating-colorblind-friendly-figures/) and saved the palette off as an Adobe Color file (`ACO`). I then took the values from the charts and mapped each party to a particular color. Then, I made ordered and proportional waffle charts using the the values and aligned colors. The results are below:

# install.packages("waffle")
# devtools::install_github("hrbrmstr/swatches")
 
library(waffle)
library(swatches)
 
national_unity <- c(`Zionist Union (27)`=27,
                    `Likud (27)`=27,
                    `Kulanu (10)`=10,
                    `Shas (7)`=7,
                    `UTJ (6)`=6)
 
right_wing <- c(`Likud (27)`=27,
                `Kulanu (10)`=10,
                `Bayit Yehudi (8)`=8,
                `Shas (7)`=7,
                `UTJ (6)`=6,
                `Yisrael Beytenu (5)`=5)
 
herzog_led <- c(`Zionist Union (27)`=27,
                `Kulanu (10)`=10,
                `Shas (7)`=7,
                `UTJ (6)`=6,
                `Meretz (5)`=5)
 
party_colors <- rev(read_aco("http://rud.is/dl/israel.aco"))
 
zion <- party_colors[1]
likud <- party_colors[2]
kulanu <- party_colors[3]
shas <- party_colors[4]
utj <- party_colors[5]
visrael <- party_colors[6]
meretz <- party_colors[7]
bayit <- party_colors[6]
 
nw <- waffle(national_unity, rows=5,
             colors=c(zion, likud, kulanu, shas, utj),
             title="\nNational unity government") +
  theme(plot.title=element_text(size=12, face="bold"))
 
rw <- waffle(right_wing, rows=5,
             colors=c(likud, kulanu, bayit, shas, utj, visrael),
             title="\nRight Wing", pad=3) +
  theme(plot.title=element_text(size=12, face="bold"))
 
hw <- waffle(herzog_led, rows=5,
             colors=c(zion, kulanu, shas, utj, meretz),
             title="\nHerzog led", pad=5) +
  theme(plot.title=element_text(size=12, face="bold"))
 
iron(nw, rw, hw)

If I knew my audience did not have color processing issues, I’d use a better palette. Regardless, these results are far better than the careless pies presented in the original story. The squares represent the same quantities in each chart and the colors also map to the parties.

Honestly, though, you could get a better idea with simple, un-tweaked base graphics bar charts:

par(mfrow=c(3,1))
barplot(national_unity, col=c(zion, likud, kulanu, shas, utj), main="National unity government")
barplot(right_wing, col=c(likud, kulanu, bayit, shas, utj, visrael), main="Right Wing")
barplot(herzog_led, col=c(zion, kulanu, shas, utj, meretz), main="Herzog led")

Please consider your readers and the message you’re trying to convey when developing visualizations, especially when you have as large an audience as the Jerusalem Post.

My disdain for pie charts is fairly well-known, but I do concede that there are times one needs to communicate parts of a whole graphically verses using just words or a table. When that need arises, I’m partial to “waffle charts” or “square pie charts”. @eagereyes did a great post a while ago on them (make sure to read the ‘debate’ between Robert and @hadleywickham in the comments, too), so head there for the low-down on them. Rather than have every waffle chart I make be a one-off creation, I made an R package for them.

There is currently one function in the package — waffle — and said function doesn’t mimic all the goodness of these charts as described in Robert’s post (yet). It does, however, do a pretty decent job covering the basics. Let’s take the oft-cited New York times “debt” graphic:

We can replicate that pretty closely in R. To make it as simple as possible, the waffle function takes a named numeric vector. If no names are specified, or you leave some names out, LETTERS will be used to fill in the gaps. The function takes your data quite literally, so if you give it a vector that sums up to, say, 10,000, then the function will try to create a ggplot object with 10,000 geom_rect elements. Needless to say, that’s a bad idea. So, I suggest using the raw numbers in the vector and passing in a scaled version of the vector to the function. That way, you can play with the values to get the desired look. Here’s the R version of of the NYT graphic:

# devtools::install_github("hrbrmstr/waffle")
library(waffle)
savings <- c(`Mortgage ($84,911)`=84911, `Auto and\ntuition loans ($14,414)`=14414, 
             `Home equity loans ($10,062)`=10062, `Credit Cards ($8,565)`=8565)
waffle(savings/392, rows=7, size=0.5, 
       colors=c("#c7d4b6", "#a3aabd", "#a0d0de", "#97b5cf"), 
       title="Average Household Savings Each Year", 
       xlab="1 square == $392")

This package evolved from a teensy gist I made earlier this year to help communicate the scope of the Anthem data breach in the US. Since then, a recent breach at Premera occurred and added to the tally. Here’s two views of that data, one with one square equalling one million people and another with one square equalling ten million people (using the blue shade from each of the company’s logos):

parts <- c(`Un-breached\nUS Population`=(318-11-79), `Premera`=11, `Anthem`=79)
 
waffle(parts, rows=8, size=1, colors=c("#969696", "#1879bf", "#009bda"), 
       title="Health records breaches as fraction of US Population", 
       xlab="One square == 1m ppl")

waffle(parts/10, rows=3, colors=c("#969696", "#1879bf", "#009bda"), 
       title="Health records breaches as fraction of US Population", 
       xlab="One square == 10m ppl"

I’m betting that gets alot bluer by the end of the year.

The function returns a ggplot object, so fonts, sizes, etc can all be customized and the source is up on github for all to play with and contribute to.

Along with adding support for filling in the chart as shown in the @eagereyes post, there will also be an htmlwidget version coming as well. Standard drill applies: issues/enhancements to github issues, feedback and your own examples in the comments.

UPDATE

Thanks to a PR by @timelyportfolio, there is now a widget option in the package.