class: inverse middle background-image: url(images/wattle_bee.jpg) background-position: 99% 99% background-size: 55% # Making inference using data plots, with application to ecological statistics ## Di Cook <br> Monash University ### vISEC <br> June 24, 2020 <br> <br> .tiny[[https://dicook.org/files/vISEC2020/slides.html](https://dicook.org/files/vISEC2020/slides.html)] <br> .footnote[Image credit: Di Cook, 2018] --- class: inverse middle center # Data plots are utilised widely in ecology, often to make decisions. They can and should be integrated into the classical statistics infrastructure. --- .large[I'm going to talk about] --
<i class="fas fa-hand-pointer fa-2x faa-float animated faa-slow " style=" color:#75A34D;"></i>
.large[.purple[inference for data plots]] --
<i class="fas fa-hand-spock fa-2x faa-wrench animated faa-slow " style=" color:#75A34D;"></i>
.large[.green[a high-throughput analysis]] --
<i class="fas fa-hand-peace fa-2x faa-vertical animated faa-slow " style=" color:#75A34D;"></i>
.large[.orange[and computer vision experiments,]] --- # 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 <img src="slides_files/figure-html/unnamed-chunk-3-1.png" width="100%" /> ] --- .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 1188 ## 2 FALSE TRUE 312 ## 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(20200624) 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(20200625) *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 a real lineup test --- # Lineup protocol I'm going to show you a page of plots -- Each has a
number
above it, this is its id -- Choose the plot that you think exhibits the
most separation
between groups -- *If you really need to choose more than one, or even not choose any, that is ok, too* --
Ready?
---
01
:
00
<br> <br> <img src="slides_files/figure-html/lineup of the wasps-1.png" width="100%" /> --- .pull-left[ The data plot is <img src="slides_files/figure-html/true wasp data plot-1.png" width="70%" /> ] .pull-right[ <br> <br> <br> <br>
My guess is that nobody picked it?
] --- .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 <img src="slides_files/figure-html/null wasp data plot-1.png" width="90%" /> ] --- class: inverse middle center Space is big, and with few data points, classes can easily be separated --
spuriously
-- <br> <br>
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 <img src="images/wasps_true.gif" style="width: 75%; align: center" /> ] .pull-right[ Top 12 PCs <img src="images/wasps_pca_true.gif" style="width: 75%; align: center" /> ] --- 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[ <img src="images/plot_turk2_300_150_5_3.png" style="width: 80%; align: center"> ] --- .left-code[ Power analysis of human evaluation relative to classical test. Effect `\(= \frac{\sqrt{n}\times |\beta|}{\sigma}\)` <br> <br> *Pooling the results from multiple people produces results that mirror the power of the classical test.* ] .right-plot[ <img src="images/majumder1.png" width="70%" /> <img src="images/majumder2.png" width="30%" /> ] --- 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
number
above it, this is its 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 <br> discrepancy <br> 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.*
<i class="fas fa-bolt faa-wrench animated faa-slow " style=" color:#75A34D;"></i>
--- class: inverse middle Two aspects of massive multiple testing - ruler on which to measure difference === .yellow[empirical Bayes] - false positives === .yellow[False Discovery Rate] -- <br>
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 --- class: inverse middle center # Bring on deep learning, and computer vision models .large[💻] --- .left-code[ .green[*Monash Masters thesis by Shuofan Zhang*] Starting from .orange[Majumder's validation study data]: `\(H_o: \beta_k = 0 ~~vs ~~H_a: \beta_k \neq 0\)` Linear vs no relationship (null) ] .right-plot[ ### Training the deep learning model Same process, but with broader range of parameter settings, and a lot more data! 200,000 samples from each of linear and null scenario generated `\(\beta_1\sim \pm U[-10, -0.1]\)` (linear, null when `\(\beta_1=0\)`) `\(\sigma \sim U[1,12]\)` `\(n=U[50,500]\)` ] --- .left-code[ ### Computer model prediction .small[ - Re-generate the 70 *data plots* using the same data in Turk study (without null plots) - Use the computer model to predict whether the 70 data plots were "linear" or "null" - The computer model's predicted accuracy over the 70 data plots are recorded as the model's performance. ] ] .right-plot[ ### Human subjects results .small[ - Calculate `\(p\)`-value associated with each lineup using the binomial formula (from Majumder), with `\(N\)`=number of evaluations and k=number of people choosing data plot - Draw conclusion: reject the null when the calculated `\(p\)`-value is smaller than `\(\alpha\)`. - The accuracy of the conclusions over the 70 lineups ] ] --- # Repeat of experiment Using same sample of `\(n\)`, `\(\beta\)`, `\(\sigma\)`, new data generated, and images created numerically by binning (to 30x30 pixels), counting and scaling counts to 0-255. Keras model fitted with 60,000 training images for each class, linear and not. Accuracy with .orange[simulated test] data, 93%. Null error 0.0179, linear error 0.1176 .footnote[Code available in the file `keras_correlation.r`] -- .pink[Its blindingly fast!] --- .left-code[ ## Accuracy <img src="images/results.png" width="100%" /> Humans beat computers. ] -- .right-plot[ ## Power analysis <img src="images/new_computer_v_human.png" width="80%" /> Humans beat computers. ] --- # Comparison of human and computer. <table> <tr> <th> </th> <th> </th> <th colspan="2"> Computer </th> </tr> <tr> <th> </th> <th> </th> <th> Not </th> <th> Linear </th> </th> </tr> <tr> <th rowspan="2"> Human </th> <th> Not </th> <td> 27 </td> <td> 0 </td> </tr> <tr> <th> Linear </th> <td> 15 </td> <td> 28 </td> </tr> </table> <br> <br> Computer tends to predict too many as "not linear". --- # 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 mange plot evaluation - Computer vision models are promising ways to scale up --- ### Additional reading .tiny[ ^ Buja et al (2009) Statistical Inference for Exploratory Data Analysis and Model Diagnostics, RSPT A <br> ^ Wickham et al (2010) Graphical Inference for Infovis, TVCG <br> ^ Hofmann et al (2012) Graphical Tests for Power Comparison of Competing Design, TVCG <br> ^ Majumder et al (2013) Validation of Visual Statistical Inference, Applied to Linear Models, JASA <br> ^ Yin et al (2013) Visual Mining Methods for RNA-Seq data: Examining Data structure, Understanding Dispersion estimation and Significance Testing, JDMGP <br> ^ Zhao, et al (2014) Mind Reading: Using An Eye-tracker To See How People Are Looking At Lineups, IJITA <br> ^ Lin et al (2015) Does host-plant diversity explain species richness in insects? Ecological Entomology<br> ^ Roy Chowdhury et al (2015) Using Visual Statistical Inference to Better Understand Random Class Separations in High Dimension, Low Sample Size Data <br> ^ Loy et al (2017) Model Choice and Diagnostics for Linear, CS <br> Mixed-Effects Models Using Statistics on Street Corners, JCGS <br> ^ Roy Chowdhury et al (2018) Measuring Lineup Difficulty By Matching Distance Metrics with Subject Choices in Crowd- Sourced Data, JCGS ] --- class: inverse middle background-image: url(images/gum.png) background-position: 99% 1% background-size: 35% # Acknowledgements Slides created via the R package [**xaringan**](https://github.com/yihui/xaringan), with **iris 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/vISEC2020/slides.html](https://dicook.org/files/vISEC2020/slides.html) and supporting files at [https://github.com/dicook/vISEC2020](https://github.com/dicook/vISEC20). <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><br />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, 2019]