---
title: Data visualisation challenge
subtitle: Harm caused by different types of drugs
author: David Selby
date: 'Wednesday 13^th^ March 2019'
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(dev.args = list(type = "cairo")) # enable anti-aliasing
```

This is an example attempt for the February 2019 Reddit
[`dataisbeautiful` data visualisation challenge][Reddit thread], to demonstrate
at [Warwick R User Group][WRUG].

## Import the data

The [original dataset][dataset] is an HTML table on Wikimedia Commons.
We can import this with the help of an RStudio add-in called **datapasta**.

```{r import}
data.frame(stringsAsFactors=FALSE,
           Drug = c("Heroin", "Cocaine",
                    "Barbiturates", "Tobacco",
                    "Alcohol", "Street methadone", "Ketamine",
                    "Amphetamine", "Benzodiazepines",
                    "4-MTA", "Buprenorphine", "Cannabis",
                    "Solvents", "Ecstasy", "Methylphenidate",
                    "GHB", "LSD", "Anabolic steroids",
                    "Khat", "Alkyl nitrites"),
           Mean.Physical.harm.plus.Mean.Dependence = c(5.78, 4.72, 4.24, 4.11, 4.08, 3.94,
                                                       3.54, 3.48, 3.46, 3.45, 3.24, 3.01,
                                                       2.91, 2.73, 2.57, 2.49, 2.36, 2.33,
                                                       1.79, 1.77),
           Mean.Physical.harm = c(2.78, 2.33, 2.23, 1.9, 2.15, 1.86,
                                  2, 1.81, 1.63, 2.15, 1.6, 1.5, 1.9,
                                  1.6, 1.32, 1.3, 1.13, 1.45, 0.75,
                                  0.93),
           Median.Physical.harm = c(2.8, 2, 2.3, 1.9, 2.15, 1.7, 2.1,
                                    1.8, 1.7, 2.15, 1.3, 1.5, 1.9, 1.6,
                                    1.3, 1.3, 1.4, 1.7, 0.75, 0.9),
           Mean.Physical.harm.NA0 = c(2.78, 2.33, 2.23, 1.24, 1.4, 1.86,
                                      2, 1.81, 1.63, 1.44, 1.6, 0.99,
                                      1.28, 1.05, 1.32, 0.86, 1.13, 1.45, 0.5,
                                      0.93),
           Acute.harm = c(2.8, 2, 2.3, 0.9, 1.9, 2.5, 2.1,
                          1.3, 1.5, 2.2, 1.2, 0.9, 2.1, 1.6,
                          1.2, 1.4, 1.7, 0.8, 0.3, 1.6),
           Chronic.harm = c(2.5, 2, 1.9, 2.9, 2.4, 1.7, 1.7,
                            1.8, 1.7, 2.1, 1.3, 2.1, 1.7, 1.6,
                            1.3, 1.2, 1.4, 2, 1.2, 0.9),
           Intravenous.harm = c(3, 3, 2.5, NA, NA, 1.4, 2.1, 2.4,
                                1.8, NA, 2.3, NA, NA, NA, 1.6, NA,
                                0.3, 1.7, NA, 0.3),
           Mean.Dependence = c(3, 2.39, 2.01, 2.21, 1.93, 2.08,
                               1.54, 1.67, 1.83, 1.3, 1.64, 1.51,
                               1.01, 1.13, 1.25, 1.19, 1.23, 0.88,
                               1.04, 0.87),
           Pleasure = c(3, 3, 2, 2.3, 2.3, 1.8, 1.9, 2,
                        1.7, 1, 2, 1.9, 1.7, 1.5, 1.4, 1.4,
                        2.2, 1.1, 1.6, 1.6),
           Psychological = c(3, 2.8, 2.2, 2.6, 1.9, 2.3, 1.7,
                             1.9, 2.1, 1.7, 1.5, 1.7, 1.2, 1.2,
                             1.3, 1.1, 1.1, 0.8, 1.2, 0.7),
           Physical = c(3, 1.3, 1.8, 1.8, 1.6, 2.3, 1, 1.1,
                        1.8, 0.8, 1.5, 0.8, 0.1, 0.7, 1,
                        1.1, 0.3, 0.8, 0.3, 0.3),
           Mean.Social.harm = c(2.54, 2.17, 2, 1.42, 2.21, 1.87,
                                1.69, 1.5, 1.65, 1.06, 1.49, 1.5,
                                1.52, 1.09, 0.97, 1.3, 1.32, 1.13,
                                0.85, 0.97),
           Intoxication = c(1.6, 1.8, 2.4, 0.8, 2.2, 1.6, 2,
                            1.4, 2, 1.2, 1.6, 1.7, 1.9, 1.2,
                            1.1, 1.4, 1.6, 1.3, 0.7, 0.8),
           Social.harm = c(3, 2.5, 1.9, 1.1, 2.4, 1.9, 1.5,
                           1.5, 1.5, 1, 1.5, 1.3, 1.5, 1, 0.8,
                           1.3, 1.3, 0.8, 1.1, 0.7),
           Healthcare.costs = c(3, 2.3, 1.7, 2.4, 2.1, 2, 1.5, 1.6,
                                1.5, 1, 1.4, 1.5, 1.2, 1.1, 1.1,
                                1.2, 1.1, 1.3, 0.8, 1.4)) -> drugs
```

## Wrangle the data

Using the R packages **dplyr** and **tidyr**---both designed for modifying and
reshaping data frames---we can filter out the variables we aren't interested in,
and assign each type of harm to one of the three broader categories of physical
harm, social harm and dependence.

```{r message = FALSE}
library(dplyr)
library(tidyr)
```

I want a lookup table that allows me to specify a type of harm and it will tell
me which of the three categories that belongs to.
Here is one way to do that, whilst avoiding repetitive code.

```{r categories}
colnames(drugs) <- gsub('\\.', ' ', colnames(drugs)) # Replace dots with spaces
categories <- setNames(rep(c('Physical', 'Dependence', 'Social'), each = 3),
                       c('Acute harm', 'Chronic harm', 'Intravenous harm',
                         'Pleasure', 'Psychological', 'Physical',
                         'Intoxication', 'Social harm', 'Healthcare costs'))
```

Now for some data wrangling.
Our first task is to get rid of the columns containing means and medians, since

a. we don't need them for now
b. we are perfectly capable of calculating them ourselves later.

The **dplyr** function `select()` and its 'helper' function `starts_with()` does
the job here, though it would also be easy to use base R functions `grep()` and
`subset()`.

From **tidyr**, `gather()` reshapes the data frame from a wide format to long,
with one row per value, rather than the original one row per drug.
The types of harm are stored as a new column called `harm`.
This structure makes it much easier to produce plots, especially in **ggplot2**
but it can help in base R as well.

```{r wrangle}
drugs_long <- drugs %>%
  as_tibble %>%
  select(-starts_with('Mean'), -starts_with('Median')) %>%
  gather('harm', 'rating', -Drug) %>%
  mutate(category = categories[harm],
         harm = factor(harm, levels = names(categories))) %>%
  na.omit()
```

The `mutate()` command performs two tasks:

1. Assign one of the three categories to each type of harm
2. Convert the column `harm` into a factor, so that we don't lose track of any
types of harm that are `NA` for some types of drug.

Let's take a peek at the new data structure.

```{r peek}
drugs_long
```

## Visualise the data

We will be using the visualisation package **ggplot2** today.

```{r message = FALSE}
library(ggplot2)
```

Many of the [entrants on Reddit][Reddit thread] are adopt *radar plots*, where
each data point is represented as the length of a spoke on a wheel.
This is straightforward to create in **ggplot2**, but first let's try something
a bit more traditional.
Here is one attempt: a stacked bar chart, facetted by category of harm.

```{r bar chart 1}
ggplot(drugs_long) +
  aes(Drug, rating, fill = harm) +
  geom_col(position = 'stack') +
  coord_flip() +
  facet_wrap(~ category)
```

It's certainly neat, but isn't so easy to draw inferences from it.
Easy improvements would be a better colour scheme, and re-ordering the drugs by
from largest to smallest total harm.
(This is left as an exercise.)

Instead, I am going for a different format, with a more modular approach.
At first glance, we can see the total harm for each drug.
A closer look reveals the relative harm by category, and finally a detailed
inspection allows the viewer to see the rating for a specific type of harm.

```{r bar chart 2, fig.width = 8, fig.height = 6}
ggplot(drugs_long) + aes(harm, rating, fill = category) +
  geom_col(position = 'stack') +
  coord_flip() +
  facet_wrap(~ Drug)
```

This is closer to what we want, but there is a lot of 'chart junk', or a low
'data to ink ratio'.
The grid lines aren't helpful, there are too many category labels, and the data
are all on a common, somewhat arbitrary scale so there is little utility in
being able to read the exact values off the \(x\)-axis.
Also, I don't like grey backgrounds amd find it annoying that `theme_grey()` is
the default colour scheme.

Let's improve things by applying **ggplot** template `theme_minimal()` as a
baseline, further customised using the `theme()` command to remove the grid
lines and \(x\)-axis labels.

```{r bar chart 3, fig.width = 8, fig.height = 6}
ggplot(drugs_long) + aes(harm, rating, fill = category) +
  geom_col(position = 'stack') +
  coord_flip() +
  facet_wrap(~ Drug) +
  labs(x = NULL, y = NULL) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text.x = element_blank())
```

One last step: launch the chart into polar coordinates and we get a radar plot!

```{r radar, fig.width = 8, fig.height = 6}
ggplot(drugs_long) + aes(harm, rating, fill = category) +
  geom_col(position = 'stack') +
  coord_polar() +
  facet_wrap(~ Drug) +
  labs(x = NULL, y = NULL) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank())
```

This doesn't tell us the categories though, so we will add a custom key.
You can combine two **ggplot2** plots into a single graphic using a function
called `gridExtra::grid.arrange()`.

```{r radar legend, fig.width = 11.5, fig.height = 7}
# Main plot
radar <- ggplot(drugs_long) + aes(harm, rating, fill = category) +
  geom_col(position = 'stack') +
  coord_polar() +
  facet_wrap(~ Drug) +
  labs(x = NULL, y = NULL) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        legend.position = 'none')

# Legend/key
key <- drugs_long %>% filter(Drug == 'Heroin') %>% mutate(rating = 1) %>%
  # ^ Pick a drug without any NAs and set all segments to the same radius
  ggplot(aes(harm, rating, label = harm)) +
  geom_col(aes(fill = category), alpha = .25) +
  # Individual segment labels
  geom_text(aes(y = rating + .3), size = 3) +
  # Broader category labels
  geom_text(aes(y = rating - .5), size = 4,
            label = as.vector(rbind('', unique(categories), ''))) +
  coord_polar() +
  labs(x = NULL, y = NULL) +
  expand_limits(y = 1.5) + # Adds outer padding
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')

# Combine into a single plot
gridExtra::grid.arrange(radar, key, ncol = 2, widths = c(.7, .3))
```

[dataset]: https://commons.wikimedia.org/wiki/File:Development_of_a_rational_scale_to_assess_the_harm_of_drugs_of_potential_misuse_(physical_harm_and_dependence,_NA_free_means).svg#Data
[WRUG]: https://www.meetup.com/Warwick-useRs/events/259364376/
[Reddit thread]: https://www.reddit.com/r/dataisbeautiful/comments/an18kt/battle_dataviz_battle_for_the_month_of_february/
