---
title: "How to make your Boxplots more informative"
output:
  flexdashboard::flex_dashboard:
    storyboard: true
    vertical_layout: fill
    theme: spacelab
    source_code: embed
    fontsize: 17pt
---

### Data Preparation

We're using data from **Chart2000.com**:
Music Charts 2000 - 2020 {data-commentary-width="500"}

```{r setup}

library(knitr)
library(kableExtra)
library(flexdashboard)
library(tidyverse)
library(ggthemes)
library(DT)
library(plotly)
library(EnvStats)
library(ggstatsplot)
library(ggtext)

knitr::opts_chunk$set(echo = FALSE)

all_songs <- read_csv(file = "chart2000-songyear-0-3-0062.csv",
                  na = c("", "-"))
attr(all_songs, "spec") <- NULL

```

```{r setup-theme}

theme_set(theme_solarized(base_size = 15))

theme_update(axis.text.x = element_text(angle = 90),
             axis.ticks.x = element_blank(),
             panel.grid.major.x = element_blank())

```

```{r data-prep}

top_artists <- all_songs %>%
  group_by(artist) %>%
  summarise(total_score = sum(indicativerevenue)) %>% 
  arrange(desc(total_score)) %>% 
  head(n = 5) %>% 
  pull(artist)

songs <- all_songs %>% 
  filter(artist %in% top_artists) %>% 
  mutate(artist = fct_infreq(artist),
         indicativerevenue = round(indicativerevenue))

datatable(songs, filter = "top")

```

------------------------------------------------------------------------

From the site [**https://chart2000.com/about.htm\#fairuse**](https://chart2000.com/about.htm#fairuse){.uri} we downloaded the file *chart2000-songyear-0-3-0062*, which you can also find on my github profile at . It contains the top 100 songs for each year from 2000 to 2020.

```{r notes-setup, ref.label = 'setup', eval = FALSE, echo = TRUE}
```

First we filter out the most successful artists: The Top 5 in terms of total score, i. e. the sum of Indicative Revenue (IR). According to , IR is *an attempt to measure the complete revenue generated by a song or album over a certain period. It does take inflation and currency conversion into account and can approximately be related to total revenue generated across the whole music chain in thousands of dollars.*

```{r notes-data-prep, ref.label = 'data-prep', eval = FALSE, echo = TRUE}
```

### First Boxplot: 
Some tweaks to a basic ggplot2 chart

using **ggthemes** by **Jeffrey B. Arnold** {data-commentary-width="600"}

```{r basic-boxplot}

ggplot(songs, aes(x = artist, y = indicativerevenue)) +
  geom_boxplot() +
  labs(x = "", y = "Indicative Revenue",
       title = "Indicative Revenue by Artist",
       subtitle = "Artists sorted by number of songs in Top 100 per year",
       caption = "Source: Chart2000.com, Songs of the year, Version 0-3-0062")

```

------------------------------------------------------------------------

This is our first attempt at summarizing indicative revenue of the Top 10 artists.

What we have already done:

-   Chosen a new theme: **ggthemes::theme_solarized()**. Thanks to package author **Jeffrey B. Arnold** and contributors
-   Increased font size
-   Ordered artists by total number of songs in top 100 for each year 2000 - 2020
-   Rotated x axis labels
-   Removed x axis tickmarks
-   Removed vertical grid lines

```{r notes-setup-theme, ref.label = 'setup-theme', eval = FALSE, echo = TRUE}
```

```{r notes-basic-boxplot, ref.label = 'basic-boxplot', eval = FALSE, echo = TRUE}
```

### 1. Show Ns Visually

by specifying **varwidth = TRUE** {data-commentary-width="600"}

```{r Ns-visually}

ggplot(songs, aes(x = artist, y = indicativerevenue)) +
  geom_boxplot(varwidth = TRUE) +
  labs(x = "", y = "Indicative Revenue",
       title = "Indicative Revenue by Artist",
       subtitle = "Artists sorted by number of songs in Top 100 per year",
       caption = "Source: Chart2000.com, Songs of the year, Version 0-3-0062")

```

------------------------------------------------------------------------

The "trick" here is simply the *varwidth = TRUE* argument in the geom_boxplot() call.

```{r notes-Ns-visually, ref.label = 'Ns-visually', eval = FALSE, echo = TRUE}
```

### 2. Show Ns Numerically 

using **EnvStats::stat_n_text()** by **Steven Millard** and **Alexander Kowarik** {data-commentary-width="600"}

```{r Ns-numerically}

ggplot(songs, aes(x = artist, y = indicativerevenue)) +
  geom_boxplot(varwidth = TRUE) +
  labs(x = "", y = "Indicative Revenue",
       title = "Indicative Revenue by Artist",
       subtitle = "Artists sorted by number of songs in Top 100 per year",
       caption = "Source: Chart2000.com, Songs of the year, Version 0-3-0062") +
  stat_n_text(y.pos = 900)

```

------------------------------------------------------------------------

We could, of course, calculate Ns and display them using *geom_text()* or *geom_label()*, but we'll use the convencience function *stat_n_text()* from the **EnvStats** package by **Steven Millard** and **Alexander Kowarik** instead.

```{r notes-Ns-numerically, ref.label = 'Ns-numerically', eval = FALSE, echo = TRUE}
```

### 3. Show Summary Statistics on Hover (Mouse-Over)

using **plotly** by **Carson Sievert** {data-commentary-width="600"}

```{r plotly}

songs %>%
  plot_ly(x = ~artist, y = ~indicativerevenue,
          type = "box")

```

***

Note that **plotly** uses a different calculation method, e. g. a different definition of outliers.

```{r notes-plotly, ref.label = 'plotly', eval = FALSE, echo = TRUE}
```


### 4. Label Outliers

via a **user-defined function** and **geom_text()** {data-commentary-width="600"}

```{r label-outliers}

is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

songs %>% 
  group_by(artist) %>% 
  mutate(outlier = ifelse(is_outlier(indicativerevenue), song, NA)) %>% 
  ggplot(aes(x = artist, y = indicativerevenue, color = artist)) +
  geom_boxplot(varwidth = TRUE) +
  geom_text(aes(label = outlier), na.rm = TRUE, nudge_y = 1500) +
  labs(x = "", y = "Indicative Revenue",
       title = "Indicative Revenue by Artist",
       subtitle = "Artists sorted by number of songs in Top 100 per year",
       caption = "Source: Chart2000.com, Songs of the year, Version 0-3-0062") +
  scale_color_discrete(guide = NULL)

```

***

A user-defined function (found on [Stackoverflow](https://stackoverflow.com/questions/33524669/labeling-outliers-of-boxplots-in-r), answer by JasonAizkalns) returns songs of outliers, NA otherwise. This newly calculated variable is passed to the **geom_text()** function.

Introducing color for a clearer group distinction. In this case we suppress the color legend, as the groups should be clear from the x axis labels.

```{r notes-label-outliers, ref.label = 'label-outliers', eval = FALSE, echo = TRUE}
```


### 5. Display Individual Data Points

using **geom_jitter()** {data-commentary-width="600"}

```{r individual-data-points}

ggplot(songs, aes(x = artist, y = indicativerevenue, color = artist)) +
  geom_boxplot(varwidth = TRUE, outlier.color = NA) +
  geom_jitter(alpha = 0.6, width = 0.2, height = 0) +
  labs(x = "", y = "Indicative Revenue",
       title = "Indicative Revenue by Artist",
       subtitle = "Artists sorted by number of songs in Top 100 per year",
       caption = "Source: Chart2000.com, Songs of the year, Version 0-3-0062") +
  scale_color_discrete(guide = NULL) +
  stat_n_text(y.pos = 900)

```

------------------------------------------------------------------------

You can get the same boxplot for different distributions, e. g. a normal distribution (most data points around the mean) vs. a u-shaped distribution (two peaks: one below and one above the mean, with very few data points near the mean).

So displaying all data points gives a better sense of the underlying distributions. It combines macro and micro levels. Inspired by **Edward Tufte**, see his great book *Envisioning Information*.

Note the use of *alpha* (opacity), *width* and *height* in **geom_jitter()**, which reduces overplotting.

```{r notes-individual-data-points, ref.label = 'individual-data-points', eval = FALSE, echo = TRUE}
```

You could use *plotly* here again to display information about songs on hover. At the moment, a hack is needed to remove outliers from the boxplot - maybe support for that will improve in a future version of *plotly*. Plotly's way seems to be to display the data points side-by-side to the boxplot, not overlay it.


### 6. Use Aesthetics to Include Information on Another Variable

Combining **shape** and **color**, using **RColorBrewer** by **Erich Neuwirth** {data-commentary-width="600"}

```{r another-variable, fig.width = 8}

songs %>% 
  rowwise() %>% 
  mutate(no1 = any(c_across(us:au) == 1, na.rm = TRUE)) %>% 
  ggplot(aes(x = artist, y = indicativerevenue)) +
    geom_boxplot(varwidth = TRUE, outlier.color = NA) +
    geom_jitter(alpha = 0.6, width = 0.2, height = 0,
                aes(shape = no1, color = no1)) +
    labs(x = "", y = "Indicative Revenue",
         title = "Indicative Revenue by Artist",
         subtitle = "Artists sorted by number of songs in Top 100 per year",
         caption = "Source: Chart2000.com, Songs of the year, Version 0-3-0062") +
    scale_color_brewer(palette = "Dark2", 
                       name = "No. 1\n(Any Country)?") +
    scale_shape_discrete(name = "No. 1\n(Any Country)?") +
    stat_n_text(y.pos = 900)

```

***

* Aesthetics for color and shape defined inside **geom_jitter()** (otherwise, there would be separate boxplots for the groups)
* Two aesthetics represented by one legend, as they refer to the same variable **and** are named the same
* Using **RColorBrewer** by **Erich Neuwirth**

```{r notes-another-variable, ref.label = 'another-variable', eval = FALSE, echo = TRUE}
```


### 7. Add Means

for comparison to medians {data-commentary-width="600"}

```{r means}

ggplot(songs, aes(x = artist, y = indicativerevenue, color = artist)) +
  geom_boxplot(varwidth = TRUE) +
  # geom_jitter(alpha = 0.6, width = 0.2, height = 0) +
  stat_summary(fun = "mean", color = "black", shape = 8) +
  labs(x = "", y = "Indicative Revenue",
       title = "Indicative Revenue by Artist",
       subtitle = "Artists sorted by number of songs in Top 100 per year",
       caption = "* Mean\n\nSource: Chart2000.com, Songs of the year, Version 0-3-0062") +
  scale_color_discrete(guide = NULL) +
  stat_n_text(y.pos = 900)

```

***

Interesting to see how in some cases (*The Black Eyed Peas*) median and mean seem to overlap, while in other cases (*Ed Sheeran*) there is a notable gap between median and mean. As is often the case, the mean is considerably higher, as it is influenced by outliers at the high end of the range of values.

Here I'd rather not show the individual data points to avoid information overflow. This chart focuses on summary statistics.

```{r notes-means, ref.label = 'means', echo = TRUE, eval = FALSE}
```


### 8. Calculate and Display Statistical Tests for Group Differences

using **ggstatsplot** by **Indrajeet Patil** {data-commentary-width="600"}


```{r ggstatsplot}

ggstatsplot::ggbetweenstats(
  data = songs,
  x = artist, xlab = "",
  y = indicativerevenue,
  ylab = "Indicative Revenue",
  plot.type = "box",
  type = "p",
  conf.level = 0.95,
  title = "Indicative Revenue by Artist"
)

```

***

The **ggstatsplot** package by **Indrajeet Patil** is very powerful and flexible, well integrated with *ggplot2* and well documented, see *help(package = "ggstatsplot")*.

Note that this package adds a large number of dependencies.

There are several tests to choose from, including non-parametric, robust and Bayesian.

See *?ggbetweenstats()* and [Website Documentation](https://indrajeetpatil.github.io/ggstatsplot/articles/web_only/ggbetweenstats.html).

```{r notes-ggstatsplot, ref.label = 'ggstatsplot', echo = TRUE, eval = FALSE}
```


### 9. Group Differences:

Example for a significant test {data-commentary-width="600"}

```{r ggstatsplot2}

songs2 <- all_songs %>% 
  filter(artist %in% c("Ed Sheeran", "Justin Timberlake", "Miley Cyrus"))

ggstatsplot::ggbetweenstats(
  data = songs2,
  x = artist, xlab = "",
  y = indicativerevenue,
  ylab = "Indicative Revenue",
  plot.type = "box",
  type = "p",
  conf.level = 0.95,
  title = "Indicative Revenue by Artist"
)

```

***

While **Ed Sheeran** reached the highest average indicative revenue, he also recorded a much higher standard deviation. Therefore, only **Justin Timberlake** averaged significantly higher than **Miley Cyrus** in this comparison. Note the p value correction for multiple comparisons (corrected Holm method).

```{r notes-ggstatsplot2, ref.label = 'ggstatsplot2', echo = TRUE, eval = FALSE}
```


### 10. Use Images as Labels

via **ggtext** by **Claus Wilke** {data-commentary-width="600"}

```{r ggtext-imagelabels, fig.width = 11, fig.height = 7}

# Download and rename files
 
labels <- c(Rihanna = "
*Rihanna*",
            Pink = "
*Pink*",
            'Maroon 5' = "
*Maroon 5*",
            'The Black Eyed Peas' = "
*Black Eyed Peas*",
            'Ed Sheeran' = "
*Ed Sheeran*"
            )

ggplot(songs, aes(x = artist, y = indicativerevenue, color = artist)) +
  geom_boxplot(varwidth = TRUE) +
  labs(x = "", y = "Indicative Revenue",
       title = "Indicative Revenue by Artist",
       subtitle = "Artists sorted by number of songs in Top 100 per year",
       caption = "Source: Chart2000.com, Songs of the year, Version 0-3-0062") +
  scale_color_discrete(guide = NULL) +
  scale_x_discrete(name = NULL, labels = labels) +
  stat_n_text(y.pos = 900) +
  theme_solarized(base_size = 14) +
  theme(axis.text.x = element_markdown(color = "black", angle = 0))

```

***

This is made possible by **Claus Wilke**'s excellent **ggtext** package.

Downloaded the files beforehand and stored them in the report / dashboard folder.

A named vector of labels is passed to *scale_x_discrete()*; besides, *theme(axis.text.x)* needs an *element_markdown()* function.

```{r notes-ggtext-imagelabels, ref.label = 'ggtext-imagelabels', echo = TRUE, eval = FALSE}
```