---
title: "Robust and Beautiful Statistical Visualization"
author: "Joses Ho and Adam Claridge-Chang"
date: "`r Sys.Date()`"
subtitle: How the estimation plot combines statistical rigour and visual design
output:
  tufte::tufte_html: default
vignette: >
  %\VignetteIndexEntry{Robust and Beautiful Statistical Visualization}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Statistical Visualization

What is data visualization? Battle-Baptiste and Rusert (2018) [^1] give a cogent and compelling definition:

> [Data visualization is] the rendering of information in a visual format to help communicate data while also generating new patterns and knowledge through the act of visualization itself.

[^1]: [W. E. B. Du Bois’s Data Portraits: Visualizing Black America](https://www.papress.com/html/product.details.dna?isbn=9781616897062). Edited by Whitney Battle-Baptiste and Britt Rusert, Princeton Architectural Press, 2018

Sadly, too many figures and visualizations in modern academic publications seemingly fail to "generate new patterns and knowledge through the act of visualization itself". Here, we propose a solution: *the estimation plot*.

# The Inadequacy of Common Plots

## The Barplot

```{r create.data1, message=FALSE, echo=FALSE, results='hide'}
library(dplyr)
set.seed(12345)

sampleN     <- 20

control_pop <- rnorm(10000, mean = 3,   sd = 0.5)
test_pop    <- rnorm(10000, mean = 3.2, sd = 0.5)

sample1     <- sample(control_pop, sampleN)
sample2     <- sample(test_pop,    sampleN)

id          <- seq(1: sampleN)
gender      <- c(rep('Male', sampleN/2), rep('Female', sampleN/2))



my.data <-
  tibble::tibble(Control = sample1, Test = sample2, 
                 ID = id, Gender = gender) %>%
  tidyr::gather(key = Group, value = Value, -ID, -Gender)


plot.frame <-
  my.data %>% 
  group_by(Group) %>% 
  summarise(mean = mean(Value),
            max  = max(Value),
            std  = sd(Value),
            sem  = sqrt(var(Value) / length(Value)),
            ci   = 1.96 * sem
            )

t.test.result <- t.test(sample1, sample2)

t.test.result$p.value # check if it is around 0.01.
```

Let's say we have performed an experiment with `r sampleN` control subjects, and `r sampleN` test subjects. We begin our data analysis by making a barplot of the data.


```{r create.themes, message=FALSE, echo=FALSE}
library(ggplot2)


estimation.plot.theme <-
  theme_classic() +
  theme(
    text                  =  element_text(family = "Work Sans"),
    axis.text             =  element_text(size = 13),
    axis.title            =  element_text(size = 15),
    
    axis.ticks.length     =  unit(7, "points"),
    axis.line.x.bottom    =  element_blank(),
    axis.ticks.x.bottom   =  element_blank(),
    axis.title.x.bottom   =  element_blank()
    )

non.floating.theme <-
  estimation.plot.theme +
  theme(
    panel.background      =  element_rect(fill = "#fffff8"),
    plot.background       =  element_rect(fill = "#fffff8"),
    legend.background     =  element_rect(fill = "#fffff8"),
    legend.box.background =  element_rect(fill = "#fffff8")
    )

floating.theme <-
  non.floating.theme +
  theme(
    axis.title.x.bottom  =  element_blank(),
    axis.ticks.x.bottom  =  element_blank()
    )

custom.fill   <- scale_fill_manual(values=c("#1F77B4", "#FF7F0E"))
custom.colour <- scale_color_manual(values=c("#1F77B4", "#FF7F0E"))

```



```{r barplot, message=FALSE, echo=FALSE, fig.width = 7, fig.height = 4}

max.U.connector    <- 5
annotate.x.pos     <- 2.7
annotate.text.size <- 5

ggplot(plot.frame, aes(x = Group, y = mean)) +
  
  estimation.plot.theme +
  
  geom_col(aes(fill = Group), width = 0.5) +
  
  scale_fill_brewer(palette="Dark2") +
  
  geom_errorbar(width = 0.,
                aes(ymin = mean - std, 
                    ymax = mean + std)) +
  
  geom_segment(x    = 1, xend = 1,
               y    = plot.frame$mean[1] + 1, 
               yend = max.U.connector) +
  geom_segment(x    = 1, xend = 2,
               y    = max.U.connector, 
               yend = max.U.connector) +
  geom_segment(x    = 2, xend = 2,
               y    = plot.frame$mean[2] + 1, 
               yend = max.U.connector) +
  
  coord_cartesian(ylim    = c(0, max.U.connector + 0.5),
                  xlim    = c(0.4, 5),
                  expand = FALSE) +
  geom_segment(x    = 0.4,  xend = 2.5,
               y    = 0,  yend = 0) +
  
  
  annotate("text", x = 1.5, y = 5.1, label = "*", vjust = 0,
           family = "Work Sans", size = 8) +
  
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 5, label = "Hides all observed values") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 4, label = "Effect size not shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 3, label = "Effect size precision not shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 2, label = "Effect size confidence interval\nnot shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 0.75, label = "Creates false dichotomy\nwith 'Significance Asterisk'") +
  
  ylab("") + guides(fill=FALSE)
  
  
```

The barplot has several shortcomings, despite enjoying widespread usage in academic journals. We're not the first ones (see [this](https://www.nature.com/articles/nmeth.2837), [this](http://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.1002128), or [that](https://onlinelibrary.wiley.com/doi/full/10.1111/ejn.13400)) to point out the myriad flaws with the barplot. Importantly, the barplot does not show us the effect size.

Alternatively, we can use a boxplot to visualize the data.

## The Boxplot

```{r boxplot, message=FALSE, echo=FALSE, fig.width = 7, fig.height = 4}
ggplot(my.data, aes(x = Group, y = Value)) +
  
  estimation.plot.theme +
  
  geom_boxplot(aes(fill = Group)) +
  
  scale_fill_brewer(palette="Dark2") +
  
  geom_segment(x    = 1, xend = 1,
               y    = plot.frame$max[1] + 0.15, 
               yend = max.U.connector) +
  geom_segment(x    = 1, xend = 2,
               y    = max.U.connector, 
               yend = max.U.connector) +
  geom_segment(x    = 2, xend = 2,
               y    = plot.frame$max[2] + 0.15, 
               yend = max.U.connector) +
  
  coord_cartesian(ylim    = c(0, max.U.connector + 0.5),
                  xlim    = c(0.4, 5),
                  expand = FALSE) +
  geom_segment(x    = 0.4,  xend = 2.5,
               y    = 0,  yend = 0) +
  
  
  annotate("text", x = 1.5, y = 5.1, label = "P < 0.05", vjust = 0,
           family = "Work Sans", size = 5) +
  
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'darkgreen', family = "Work Sans",
           y = 5, label = "Medians, quartiles,") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'darkgreen', family = "Work Sans",
           y = 4.7, label = "minima, and maxima shown,") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 4.4, label = "but all observations not shown") +
  
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 3.75, label = "Effect size not shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 3, label = "Effect size precision not shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 2, label = "Effect size confidence interval\nnot shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 1, label = "Creates false dichotomy\nwith P value threshold") +
  
  ylab("") + guides(fill=FALSE)
```

Unfortunately, the boxplot still doesn't show all our data. We still lack information about the underlying distribution of your data. Is it normally distributed? Is there skew in the points? What is the sample size? More importantly, boxplots do not display the effect size.

To display several data points across one or more categories, we can use the jitter plot.

## The Jitter Plot

```{r jitter.plot, message=FALSE, echo=FALSE, fig.width = 7, fig.height = 4}

pval.round <- signif(t.test.result$p.value, 2)

ggplot(my.data, aes(x = Group, y = Value)) +
  
  estimation.plot.theme +
  
  geom_jitter(aes(colour = Group), width = 0.25) +
  
  scale_color_brewer(palette="Dark2") +
  
  geom_segment(x    = 1, xend = 1,
               y    = plot.frame$max[1] + 0.15, 
               yend = max.U.connector) +
  geom_segment(x    = 1, xend = 2,
               y    = max.U.connector, 
               yend = max.U.connector) +
  geom_segment(x    = 2, xend = 2,
               y    = plot.frame$max[2] + 0.15, 
               yend = max.U.connector) +
  
  coord_cartesian(ylim    = c(0, max.U.connector + 0.5),
                  xlim    = c(0.4, 5),
                  expand = FALSE) +
  geom_segment(x    = 0.4,  xend = 2.5,
               y    = 0,  yend = 0) +
  
  
  annotate("text", x = 1.5, y = 5.1, size = 4, vjust = 0,
           family = "Work Sans",
           label = stringr::str_interp("P = ${pval.round}")) +
  
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'darkgreen', family = "Work Sans",
           y = 5, label = "All observed values shown,") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 4.7, label = "but underlying distribution,") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 4.4, label = "not accurately depicted") +
  
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 3.75, label = "Effect size not shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 3, label = "Effect size precision not shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 2, label = "Effect size confidence interval\nnot shown") +
  annotate("text", x = annotate.x.pos, size = annotate.text.size, 
           hjust = 0, color = 'red', family = "Work Sans",
           y = 1, label = "Creates false dichotomy\nwith P value threshold") +
  
  ylab("") + guides(colour=FALSE)
```

Jitter plots avoid overlapping datapoints (i.e. datapoints with the same y-value) by adding a random factor to each point along the orthogonal x-axes. Thus, while a jitter plot displays all datapoints (implicitly indicating the sample size visually), it might not accurately depict the underlying distribution of the data.

# Introducing the Estimation Plot 

```{r estimation.plot, warning=FALSE, message=FALSE, echo=FALSE, fig.width = 6, fig.height = 8}
library(dabestr)


gardner.altman <- 
  dabest(my.data, x = Group, y = Value, 
         idx = c("Control", "Test"),
         paired = FALSE) %>%
  plot(palette = "Dark2",
       theme = estimation.plot.theme,
       rawplot.type = 'swarmplot',
       width = 0.1
       # rawplot.markersize = 6,
       # effsize.markersize = 8
       )


annotation <- 
  ggplot() +
  estimation.plot.theme +
  theme(axis.text  = element_blank(),
        axis.line  = element_blank(),
        axis.ticks = element_blank()) +
  coord_cartesian(ylim   = c(0.5, 5.5),
                  xlim   = c(0, 2),
                  expand = FALSE) +
  
  annotate("text", x = 0.1, size = annotate.text.size, 
           hjust = 0, color = 'darkgreen', family = "Work Sans",
           y = 5, label = "All observed values shown")+

  annotate("text", x = 0.1, size = annotate.text.size, 
           hjust = 0, color = 'darkgreen', family = "Work Sans",
           y = 4, label = "Effect size is shown") +
  
  annotate("text", x = 0.1, size = annotate.text.size, 
           hjust = 0, color = 'darkgreen', family = "Work Sans",
           y = 3, label = "Effect size precision is displayed.") +
  
  annotate("text", x = 0.1, size = annotate.text.size, 
           hjust = 0, color = 'darkgreen', family = "Work Sans",
           y = 2, label = "Confidence and likelihood of effect size are shown") +
  
  annotate("text", x = 0.1, size = annotate.text.size, 
           hjust = 0, color = 'darkgreen', family = "Work Sans",
           y = 1, label = "No significance testing shown, so no false dichotomy") +
  
  ylab("") + guides(colour=FALSE)

cowplot::plot_grid(gardner.altman, annotation, nrow = 2)

```

Shown above is a [Gardner-Altman estimation plot](https://en.wikipedia.org/wiki/Estimation_statistics#Gardner-Altman_plot)[^2]. This plot has two key features. Firstly, it presents all datapoints as a swarmplot, which orders each point to display the underlying distribution. Secondly, an estimation plot presents the effect size as a _bootstrap 95% confidence interval_ (95% CI) on a _separate but aligned axes_. where the effect size is displayed to the right of the war data, and the mean of the test group is aligned with the effect size. *Thus, estimation plots are robust, beautiful, and convey important statistical information elegantly and efficiently*.


[^2]: The plot draws its name from [Martin J. Gardner](https://www.independent.co.uk/news/people/obituary-professor-martin-gardner-1470261.html) and [Douglas Altman](https://www.bmj.com/content/361/bmj.k2588), who are credited with [creating the design](https://www.bmj.com/content/bmj/292/6522/746.full.pdf) in 1986.


An estimation plot obtains and displays the 95% CI through *nonparametric bootstrap resampling*. This enables visualization of the confidence interval as a graded sampling distribution. This is different from the original Gardner-Altman design: the 95% CI was computed through parametric methods, and displayed as a vertical error bar.

<!-- You can read more about bootstrap resampling in this [vignette](bootstrap-confidence-intervals.html).  -->

## Estimation Statistics

Estimation plots emerge from **estimation statistics** a simple [framework](https://thenewstatistics.com/itns/) that avoids the [pitfalls of significance testing](https://www.nature.com/articles/nmeth.3288). It focuses on the effect sizes of one's experiment/interventions, and uses familiar statistical concepts: means, mean differences, and error bars.

Significance testing calculates the probability (the *P* value) that the experimental data would be observed, if the intervention did not produce a change in the metric measured (i.e. the null hypothesis). This leads analysts to apply a false dichotomy on the experimental intervention.

Estimation statistics, on the other hand, focuses on the magnitude of the effect (the effect size) and its precision. This encourages analysts to gain a deeper understanding of the metrics used, and how they relate to the natural processes being studied.

## An Estimation Plot For Every Experimental Design

For each of the most routine significance tests, **there is an estimation replacement**:

*Unpaired Student’s t-test → Two-group estimation plot*

```{r gardner.altman.unpaired, warning=FALSE, message=FALSE, echo=FALSE, fig.width = 6, fig.height = 4}


two.group.unpaired <- 
  my.data %>%
  dabest(Group, Value, 
         # The idx below passes "Control" as the control group, 
         # and "Group1" as the test group. The mean difference
         # will be computed as mean(Group1) - mean(Control1).
         idx = c("Control", "Test"), 
         paired = FALSE)

plot(two.group.unpaired, theme = estimation.plot.theme, color.column = Gender)

``` 

*Paired Student’s t-test → Paired estimation plot*

The Gardner-Altman estimation plot can also display effect sizes for repeated measures (*aka* a paired experimental design) using a [Tufte slopegraph](http://charliepark.org/slopegraphs/) instead of a swarmplot.

```{r gardner.altman.paired, warning=FALSE, message=FALSE, echo=FALSE, fig.width = 6, fig.height = 4}
two.group.paired <- 
  my.data %>%
  dabest(Group, Value, 
         idx = c("Control", "Test"), 
         paired = TRUE, id.col = ID)

plot(two.group.paired, theme = estimation.plot.theme, color.column = Gender)
``` 


```{r create.data2, message=FALSE, echo=FALSE}
set.seed(12345)

N = 20
control1_pop <- rnorm(10000, mean = 3,   sd = 0.5)
control2_pop <- rnorm(10000, mean = 3,   sd = 1  )

test1_pop    <- rnorm(10000, mean = 3.2, sd = 0.5)
test2_pop    <- rnorm(10000, mean = 2.8, sd = 0.5)
test3_pop    <- rnorm(10000, mean = 3.2, sd = 1.5)
test4_pop    <- rnorm(10000, mean = 2.8, sd = 1  )
test5_pop    <- rnorm(10000, mean = 2.8, sd = 0.5)

c1     <- sample(control1_pop, sampleN)
c2     <- sample(control2_pop, sampleN)
g1     <- sample(test1_pop,    sampleN)
g2     <- sample(test2_pop,    sampleN)
g3     <- sample(test3_pop,    sampleN)
g4     <- sample(test4_pop,    sampleN)
g5     <- sample(test5_pop,    sampleN)


id          <- 1: sampleN
gender      <- c(rep('Male', sampleN/2), rep('Female', sampleN/2))


wide.data <- 
  tibble::tibble(
    Control1 = c1, Control2 = c2,
    Group1 = g1, Group2 = g2, Group3 = g3, Group4 = g4, Group5 = g5,
    Gender = gender, ID = id)


tidy.data   <- 
  wide.data %>%
  tidyr::gather(key = Group, value = Value, -ID, -Gender)
```

For comparisons between 3 or more groups that typically employ analysis of variance (ANOVA) methods, one can use the [Cumming estimation plot](https://en.wikipedia.org/wiki/Estimation_statistics#Cumming_plot)[^3], which can be considered a variant of the Gardner-Altman plot. 

[^3]: The Cumming plot is named after [Geoff Cumming](https://www.youtube.com/watch?v=nDN-hcKR7j8), and draws its design heavily from his 2012 textbook [Understanding the New Statistics](https://www.routledge.com/Understanding-The-New-Statistics-Effect-Sizes-Confidence-Intervals-and/Cumming/p/book/9780415879682). 

*One-way ANOVA + multiple comparisons → Multi two-group estimation plot*

```{r multi.two.group.unpaired, warning=FALSE, message=FALSE, echo=FALSE, fig.width = 9, fig.height = 5}

multi.two.group.unpaired <- 
  tidy.data %>%
  dabest(Group, Value, 
         idx = list(c("Control1", "Group1"), 
                    c("Control2", "Group2"),
                    c("Group3", "Group4")),
         paired = FALSE
  )

plot(multi.two.group.unpaired, theme = estimation.plot.theme,
     tick.fontsize = 13, axes.title.fontsize = 15,
     color.column = Gender)
```

The effect size and 95% CIs are still plotted a separate axes, but unlike the Gardner-Altman plot, this axes is positioned beneath the raw data. Such a design frees up visual space in the upper panel, allowing the display of  summary measurements (mean ± standard deviation) for each group. These are shown as gapped lines to the right of each group. The mean of each group is indicated as a gap in the line, adhering to Edward Tufte's dictum to [keep the data-ink ratio low](https://medium.com/@plotlygraphs/maximizing-the-data-ink-ratio-in-dashboards-and-slide-deck-7887f7c1fab).

*Repeated measures ANOVA → Multi paired estimation plot*

```{r multi.two.group.paired, warning=FALSE, message=FALSE, echo=FALSE, fig.width = 9, fig.height = 5}

multi.two.group.paired <- 
  tidy.data %>%
  dabest(Group, Value, 
         idx = list(c("Control1", "Group1"), 
                    c("Control2", "Group2"),
                    c("Group3", "Group4")),
         paired = TRUE, id.col = ID
  )

plot(multi.two.group.paired, theme = estimation.plot.theme,
     tick.fontsize = 13, axes.title.fontsize = 15,
     color.column = Gender)
```

Ordered groups ANOVA → Shared-control estimation plot

```{r shared.control, warning=FALSE, message=FALSE, echo=FALSE, fig.width = 9, fig.height = 5}

shared.control <- 
  tidy.data %>%
  dabest(Group, Value, 
         idx = c("Control1", "Group1", "Group2", "Group3", "Group4"),
         paired = FALSE
  )

plot(shared.control, theme = estimation.plot.theme, 
     tick.fontsize = 13, axes.title.fontsize = 15,
     color.column = Gender)

```

# Estimation Plots: The Way Forward

In summary, estimation plots offer five key benefits relative to conventional plots:

|                                 | Barplot  | Boxplot | Jitter plot | Estimation Plot |
|:--------------------------------|----------|---------|-------------|-----------------|
| Avoid false dichotomy           | ✘        | ✘       | ✘           | ✔              | 
| Display all observed values     | ✘        | ✘       | ✘           | ✔              | 
| Focus on effect size            | ✘        | ✘       | ✘           | ✔              | 
| Visualize effect size precision | ✘        | ✘       | ✘           | ✔              | 
|Show mean difference distribution| ✘        | ✘       | ✘           | ✔              | 

You can create estimation plots using the DABEST (Data Analysis with Bootstrap Estimation) packages, which are available in [Matlab](https://github.com/ACCLAB/DABEST-Matlab),  [Python](https://github.com/ACCLAB/DABEST-python), and [R](https://github.com/ACCLAB/dabestr).

