Statistics on Street Corners: Conducting Inference for Data Plots

class: inverse middle
background-image: url(images/wattle_bee.jpg)
background-position: 99% 99%
background-size: 55%

# Statistics on Street Corners: Conducting Inference for Data Plots

## Di Cook Monash University

### NCB2021 June 22, 2021

.tiny[[https://dicook.org/files/NCB2021/slides.html](https://dicook.org/files/NCB2021/slides.html)]
 
.footnote[Image credit: Di Cook, 2018]

---
class: middle center inverse

# Hello 👋🏻
--

<center><img src="http://dicook.org/img/dicook-2019.png" style="width:25%; border-radius: 50%;"> </center>
--

Professor, Monash University, Melbourne, Australia

---
class: inverse middle center

# Data plots are often to make decisions, especially for deciding on whether a model fits. They can and should be integrated into the classical statistics infrastructure.

---

.large[I'm going to talk about]

.large[.purple[inference for data plots]]

.large[.green[examples, deer sightings in the wild, gene expression, making disease maps, and RNA-seq]]

.large[.orange[and how you can do this too.]]

---
# Inference for data plots requires

1. the plot is a statistic 
--

2. the type of plot (specified by a grammar) implicitly defines the null hypothesis
--

3. a null generating mechanism provides draws from the sampling distribution, among which to embed the data plot
--

4. human (computer) observers are engaged to conduct a lineup test
--

5. statistical significance and power can be computed based on the proportion of observers choosing the data plot from the lineup

---
# Why is a plot a statistic?

Many of you (hopefully) use `ggplot2` to make your plots with a grammar of graphics.

```r
ggplot(data=DATA) + 
  geom_something(
*   mapping=aes(x=VAR1, y=VAR2, colour=VAR3)
  ) +
  extra nice styling
```
--

*A statistic is a function of a random variable(s).*  This is how the mapping can be interpreted.

---

.left-code[

Adding data gives a visual statistic

```r
# Get some data
library(amt)
data("deer")
data("sh_forest")
rsf1 <- deer %>% 
 random_points(n=1500) %>% 
 extract_covariates(sh_forest) %>% 
 mutate(forest = sh.forest == 1) %>%
 rename(x=x_, y=y_, sighted=case_)

# Plot it
*ggplot(data=rsf1) +
  geom_point(
*   aes(x=x, y=y, colour=sighted),
    alpha=0.7) +
  extra nice styling
```

]
.right-plot[

Observed value of the statistic

]

---

.left-full[

```r
*ggplot(rsf1) +
  geom_bar(
*   aes(x=sighted, fill=forest),
    position = "fill") + 
  extra nice styling
```

For sighted vs forest habitat the *mapping requires call to* `stat=count`:

```
## # A tibble: 4 x 3
## # Groups: sighted [2]
## sighted forest count
## <lgl> <lgl> <int>
## 1 FALSE FALSE 1208
## 2 FALSE TRUE 292
## 3 TRUE FALSE 560
## 4 TRUE TRUE 266
```

]

.right-plot[

Observed value of statistic

<img src="slides_files/figure-html/unnamed-chunk-6-1.png" width="100%" />
]

---
# Null generating mechanism: Example 1

.left-code[

What's the null? What would be uninteresting?

```r
ggplot(DATA) + 
  geom_POINT(
*   aes(x=x, y=y, colour=sighted),
    alpha=0.7) +
  extra nice styling
```
]
--

.right-plot[

`\(H_o:\)` Sightings are uniformly distributed in space

`\(H_a:\)` Sightings are NOT uniformly distributed in space

Null generating mechanism could be to permute the labels of sighted variable. (Or could simulated a second uniform set of points.)
]

---
# Null generating mechanism: Example 2

.left-code[

What's the null? What would be uninteresting?

```r
ggplot(DATA) + 
  geom_BAR(
*   aes(x=sighted, fill=forest),
    position = "fill") + 
  extra nice styling
```
]
--

.right-plot[

`\(H_o:\)` No relationship between sighted and forest habitat

`\(H_a:\)` Sightings in forest habitat more likely

Null generating mechanism could also be permute the labels of sighted (or forest) variable. (Or could simulate from a binomial.)
]
---
class: inverse middle center

# Pretend you haven't seen the data plot

---

.left-code[

Which plot is different from the rest?

```r
set.seed(20210622)
library(nullabor)
*l <- lineup(null_permute("sighted"),
* rsf1, n=6)
ggplot(l) + 
 geom_point(
 aes(x=x, y=y, colour=sighted), 
 alpha=0.3) +
 facet_wrap(~.sample, ncol=2) + 
 extra nice styling 
```

]

.right-plot[
<img src="slides_files/figure-html/unnamed-chunk-10-1.png" width="100%" />

]

You say 1? Oh, that is the data plot.

---

.left-code[

```r
set.seed(20210622)
*l <- lineup(null_permute("sighted"),
* rsf1, n=9)
ggplot(l) +
 geom_bar(
 aes(x=sighted, fill=forest), 
 position = "fill") + 
 facet_wrap(~.sample, ncol=3) + 
 extra nice styling
```

*In which plot is the light brown bar on the right the tallest?*

]
.right-plot[
<img src="slides_files/figure-html/unnamed-chunk-12-1.png" width="100%" />

]

Did you say 5? You're good!

---
class: middle center

# In each case, the data plot was identifiable, and the null hypothesis would be rejected

---
class: inverse middle center

# Inference for graphics infrastructure

---
background-image: url(images/vis_inf.png)
background-size: contain

---
class: inverse middle center

# Visual inference broadens the scope of statistics

---
class: inverse middle center

# Let's do an actual lineup test

---
# Lineup protocol

I'm going to show you a page of plots

Each has a .orange[number] above it, this is its .orange[id]

Choose the plot that you think exhibits the .orange[most separation] between groups

*If you really need to choose more than one, or even not choose any, that is ok, too*

## Ready?

---

---

.pull-left[
The data plot is

<img src="slides_files/figure-html/true wasp data plot-1.png" width="70%" />
]
.pull-right[

## My guess is that you didn't picked this one?

]

---

.pull-left[
<img src="images/toth_PRSB_2010.png" style="width: 80%; align: center" />

> LDA resulted in ... that gynes had the most divergent expression patterns

.footnote[Toth et al (2010) Proc. of the Royal Society]
]

.pull-right[
<img src="images/toth_Science_2007.png" style="width: 75%; align: center" />

> ... show that foundress and worker brain profiles are more similar to each other than to the other groups.

.footnote[Toth et al (2007) Science]
]
---

.pull-left[
True data

<img src="slides_files/figure-html/true wasp data plot again-1.png" width="90%" />
]
.pull-right[

Null data

]
---
class: inverse middle center

Space is big, and with few data points, classes can easily be separated

.orange[spuriously].

## The lineup protocol can help people understand the problem.

---

If you first do dimension reduction (e.g. PCA), and then LDA, the problem goes away. LDA into three dimensions shown below.

.pull-left[

All data

]

.pull-right[

Top 12 PCs

]

---
class: inverse middle center

# Let's do another actual lineup test

---
# Lineup protocol

I'm going to show you a page of plots

Each has a .orange[number] above it, this is its .orange[id]

Choose the plot that you think exhibits the .orange[most separation] between groups

*If you really need to choose more than one, or even not choose any, that is ok, too*

## Ready?

---
background-color: white

---
class: inverse middle center

# Try another one

---
background-color: white

---
# Answers

- In both of these lineups the data plot is in position 8. 
- The first display is a .orange[choropleth map], where big spatial areas dominate.
- The second display is a .orange[hexagon tile map], a new display that we have designed for Australian disease communication.

---
background-color: white
# Thyroid cancer incidence

.pull-left[
<center>
<img src="images/choro-1.png" style="width: 100%; align: center" /> </center>

]
.pull-right[

]

The work is motivated by the Australian Cancer Atlas. Our experiment, using lineups showed that the hexagon tile map more effectively communicates the distribution of cancer incidence.

.tiny[Kobakian and Cook (unpublished) https://github.com/srkobakian/experiment]
---
class: inverse middle center

# What's that you say? That people can't look at so many plots?

Crowd-sourcing can help here.

---
# Validation experiment

[Majumder et al (2013)](https://www.tandfonline.com/doi/abs/10.1080/01621459.2013.808157) conducted validation study to compare the performance of the lineup protocol, assessed by human evaluators, in comparison to the classical test, using subjects employed with Amazon's Mechanical Turk.

---
## Explanation of experiment

Read about it at http://datascience.unomaha.edu/turk/exp2/index.html

.left-code[

`\(H_o: \beta_k = 0 ~~vs ~~H_a: \beta_k \neq 0\)`

- 70 lineups of size 20 plots: 
    - `\(n=100, 300\)` 
    - `\(\beta \in [-6, 4.5]\)` 
    - `\(\sigma=5, 12\)`
- 351 evaluations by human subjects

]

.right-plot[

]

---

.left-code[

Power analysis of human evaluation  relative to classical test.

Effect `\(= \frac{\sqrt{n}\times |\beta|}{\sigma}\)`

*Pooling the results from multiple people produces results that mirror the power of the classical test.*

]
.right-plot[

]
---
class: inverse middle center

# High-throughput analysis

.large[😓]

The wasps example made us worried about our own RNA-Seq analyses!

---
# Lineup of our own data

I'm going to show you a page of plots

Each has a .orange[number] above it, this is its orange[id].

Choose the plot that you think exhibits the 
- steepest green line
- with relatively small spread of the green points

## Ready?

---
background-image: url(images/plot_turk9_interaction_2_1.png)
background-position: 50% 85%
background-size: 75%

---
background-image: url(images/RNASeq_disagreement.png)
background-position: 90% 15%
background-size: 40%

Experimental design 2x2 factorial:

- Two genotypes (EV, RPA)
- Two growing conditions (I, S)
- Three reps for each treatment
- Approx 60,000 genes

Results from two different procedures, edgeR and DESeq provided conflicting numbers of significant genes, but on the order of 300 significant genes.

One of the top genes was selected for the lineup study, and independent observers engaged through Amazon's Mechanical Turk. 
---
background-image: url(images/RNASeq_explanation.png)
background-position: 90% 60%
background-size: 60%

## How does a discrepancy happen?

---
# Turk results

Is there any significant structure in our data?

- 24 lineups were made, only one shown to an observer
- 5 different positions of the data plot
- 5 different sets of null plots

Pooling results gave a detection rate of 0.65, which is high. *There is some structure to our data.*

---
class: inverse middle

Two aspects of massive multiple testing

- ruler on which to measure difference === .yellow[empirical Bayes]
- false positives === .yellow[False Discovery Rate]

Even with these, mistakes can happen, and visualising the data remains valuable.

---

background-image: url(images/RNASeq_top25.png)
background-position: 50% 50%
background-size: contain

---
background-image: url(http://dicook.github.io/nullabor/reference/figures/nullabor_hex.png)
background-position: 70% 1%
background-size: 15%

# How to do this yourself

Get a copy of the `nullabor` package

```r
install.packages("nullabor")
```

```r
# install.packages("remotes")
remotes::install_github("dicook/nullabor")
```

Look at the "Get started" documentation at http://dicook.github.io/nullabor/index.html

---
# Thanks for listening!

Here's what I hope you heard:

- Plots can be embedded into an inferential framework
- This extends the applicability of statistics to more complex problems
- Crowd-sourcing can help manage plot evaluation

---
### Additional reading

.tiny[
^ Buja et al (2009) Statistical Inference for Exploratory Data Analysis and Model Diagnostics, RSPT A 
^ Wickham et al (2010) Graphical Inference for Infovis, TVCG 
^ Hofmann et al (2012) Graphical Tests for Power Comparison of Competing Design, TVCG 
^ Majumder et al (2013) Validation of Visual Statistical Inference, Applied to Linear Models, JASA 
^ Yin et al (2013) Visual Mining Methods for RNA-Seq data: Examining Data structure, Understanding Dispersion estimation and Significance Testing, JDMGP 
^ Zhao, et al (2014) Mind Reading: Using An Eye-tracker To See How People Are Looking At Lineups, IJITA 
^ Lin et al (2015) Does host-plant diversity explain species richness in insects? Ecological Entomology 
^ Roy Chowdhury et al (2015) Using Visual Statistical Inference to Better Understand Random Class Separations in High Dimension, Low Sample Size Data, CS 
^ Loy et al (2017) Model Choice and Diagnostics for Linear, Mixed-Effects Models Using Statistics on Street Corners, JCGS 
^ Roy Chowdhury et al (2018) Measuring Lineup Difficulty By Matching Distance Metrics with Subject Choices in Crowd- Sourced Data, JCGS 
^ Vanderplas et al (2020) Testing Statistical Charts: What Makes a Good Graph? ARSIA 
^ Vanderplas et al (2021) Statistical significance calculations for scenarios in visual inference. Stat. 
]

---
class: inverse middle
background-image: url(images/wattle_bee.jpg)
background-position: 99% 1%
background-size: 25%

# Acknowledgements

Slides created via the R package [**xaringan**](https://github.com/yihui/xaringan), with **wattle theme** created from [xaringanthemer](https://github.com/gadenbuie/xaringanthemer).

The chakra comes from [remark.js](https://remarkjs.com), [**knitr**](http://yihui.name/knitr), and [R Markdown](https://rmarkdown.rstudio.com).

Slides are available at [https://dicook.org/files/NCB2021/slides.html](https://dicook.org/files/NCB2021/slides.html) and supporting files at [https://github.com/dicook/NCB2021](https://github.com/dicook/NCB2021).

<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png" /></a> This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a>.

.footnote[Image credit: Di Cook, 2018]