Visual Inference

.info-box.w-50.bg-white[
These slides are viewed best by Chrome or Firefox and occasionally need to be refreshed if elements did not load properly. See <a href="visinf-slides.pdf">here for the PDF </a>. 
]

---

background-size: cover
class: title-slide
count: false

# .monash-blue[Visual Inference]

<h2 style="font-weight:900!important;">絵による推論</h2>

.bottom_abs.width100[

Presenter: *Emi Tanaka*

Department of Econometrics and Business Statistics, &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Monash University, Melbourne, Australia

emi.tanaka@monash.edu

<a href="https://twitter.com/statsgen" style="color:black"> @statsgen</a>

27 Mar 2021 @ Fukuoka R

]

---

.footnote.monash-black.f4[
 スライドは <a href="https://emitanaka.org/slides/FukuokaR2021" style="color:#6F7C4D!important;">emitanaka.org/slides/FukuokaR2021</a>
]

---

# A little about me

# 自己紹介

* Monash大学で統計学の助教授をやってます
--

.flex[
.w-50[
<div id="monash" style="width:432px;height:288px;" class="leaflet html-widget"></div>
<script type="application/json" data-for="monash">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"setView":[[-37.9126,145.1327],4,[]],"calls":[{"method":"addTiles","args":["//{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",null,null,{"minZoom":0,"maxZoom":18,"tileSize":256,"subdomains":"abc","errorTileUrl":"","tms":false,"noWrap":false,"zoomOffset":0,"zoomReverse":false,"opacity":1,"zIndex":1,"detectRetina":false,"attribution":"© <a href=\"http://openstreetmap.org\">OpenStreetMap<\/a> contributors, <a href=\"http://creativecommons.org/licenses/by-sa/2.0/\">CC-BY-SA<\/a>"}]},{"method":"addMarkers","args":[-37.9126,145.1327,null,null,null,{"interactive":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},"Monash University",null,null,null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]}],"limits":{"lat":[-37.9126,-37.9126],"lng":[145.1327,145.1327]}},"evals":[],"jsHooks":[]}</script>

]
{{content}}

]

<div class="w-50">
<div id="japan" style="width:432px;height:288px;" class="leaflet html-widget"></div>
<script type="application/json" data-for="japan">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"calls":[{"method":"addTiles","args":["//{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",null,null,{"minZoom":0,"maxZoom":18,"tileSize":256,"subdomains":"abc","errorTileUrl":"","tms":false,"noWrap":false,"zoomOffset":0,"zoomReverse":false,"opacity":1,"zIndex":1,"detectRetina":false,"attribution":"© <a href=\"http://openstreetmap.org\">OpenStreetMap<\/a> contributors, <a href=\"http://creativecommons.org/licenses/by-sa/2.0/\">CC-BY-SA<\/a>"}]},{"method":"addCircleMarkers","args":[36.6485,138.195,10,null,null,{"interactive":true,"className":"","stroke":true,"color":"#03F","weight":5,"opacity":0.5,"fill":true,"fillColor":"#03F","fillOpacity":0.2},null,null,"長野",null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]},{"method":"addCircleMarkers","args":[33.2631,130.3009,10,null,null,{"interactive":true,"className":"","stroke":true,"color":"#03F","weight":5,"opacity":0.5,"fill":true,"fillColor":"#03F","fillOpacity":0.2},null,null,"佐賀",null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]},{"method":"addMarkers","args":[33.5902,130.4017,null,null,null,{"interactive":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},"福岡",null,null,null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]}],"limits":{"lat":[33.2631,36.6485],"lng":[130.3009,138.195]}},"evals":[],"jsHooks":[]}</script>

</div>

* 父親は佐賀市、母親は長野県からです
--

* 統計学とRを教える気歴は~10年です（±３年standard deviation)

---

# .yellow[Visual Inference]

## .yellow[絵による推論]

.tl.ml7.f1[

…の前に、

<ul class="fa-ul">
 <li>hypothesis testing (frequentist version)</li>
 <li>linear model</li>
 <li>residual analysis</li>
 <li>null plots{{content}}</li>

</ul>
]

with !

---

# レッスン .circle-big[一]

.tl.ml7.f1[

<ul class="fa-ul">
 <li>hypothesis testing (frequentist version)</li>
 <li style="opacity:0.3;">linear model</li>
 <li style="opacity:0.3;">residual analysis</li>
 <li style="opacity:0.3;">null plots{{content}}</li>

</ul>
]

---

# Testing coin bias

* 私はこの１００円持ってます　<img src="images/coin-tails.jpg" height = "50px" style="vertical-align:middle;">
* 私はこのコインは表に偏ってると思ってます。
* 即ち、５０％以上の確率で表がでるのでしょうか？

## 実験その１

コインを１０回フリップします。

結果：

７回は表、３回は裏でした。

コインは表に偏ってますか？

---

## 実験その２

コインを１００回フリップします。

結果：

７０回は表、３０回は裏でした。また７割は表でした。

コインは表に偏ってますか？

---

# (Frequentist) hypotheses testing framework

* Suppose `$X$` is the number of heads out of `$n$` independent tosses.
* Let `$p$` be the probability of getting a head for this coin.

| | |
|-----|-----|
|**Hypotheses** | `$H_0: p = 0.5$` vs. `$H_1: p > 0.5$` |
|**Assumptions** | Each toss is independent with equal chance of getting a head. |
|**Test statistic** | `$X \sim B(n, p)$` under `$H_0$`. |
|**P-value** .font_small[(or critical value or confidence interval)] | `$P(X \geq x)$` where `$x$` is the observed test statistic. |
|** Conclusion** | Reject null hypothesis when the `$p$`-value is less than some significance level `$\alpha$`. Usually `$\alpha = 0.05$`.|

* 実験その１の p-value は `$P(X \geq 7) \approx 0.17$` .inline-code[1 - pbinom(6, 10, 0.5)] Fail to reject `$H_0$`.  
* 実験その２の p-value は  `$P(X \geq 70) \approx 0.00004$` .inline-code[1 - pbinom(69, 100, 0.5)] Reject `$H_0$`.

<a href="#comparison" class="button">Inference framework</a>

---

# レッスン .circle-big[二]

.tl.ml7.f1[

<ul class="fa-ul">
 <li style="opacity:0.3;">hypothesis testing (frequentist version)</li>
 <li>linear model</li>
 <li style="opacity:0.3;">residual analysis</li>
 <li style="opacity:0.3;">null plots{{content}}</li>

</ul>
]

---

# Black Cherry Trees: Volume vs Girth

```r
library(ggplot2)
ggplot(trees, aes(log(Girth), log(Volume))) +
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
```

]
.w-60[

.monash-bg-blue.white.border-box.w-100[
Consider the model:

`$$\log(\texttt{Volume}_i) = \beta_0 + \beta_1 \log(\texttt{Girth}_i) + e_i$$`
where `$e_i \sim NID(0, \sigma^2)$`.

]

Note: *NID* means normal, independent and identically distributed.

]

---

# Model fit

.h-80.overflow-scroll[

```r
library(magrittr) # for `%>%`
library(broom) # for `augment`
fit <- lm(log(Volume) ~ log(Girth), data = trees)
res <- augment(fit)
tidy(fit)
```

```
## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2.35 0.231 -10.2 4.18e-11
## 2 log(Girth) 2.20 0.0898 24.5 6.36e-21
```

```r
res
```

```
## # A tibble: 31 x 7
## `log(Volume)` `log(Girth)` .fitted .std.resid .hat .sigma .cooksd
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2.33 2.12 2.30 0.281 0.151 0.117 0.00703 
## 2 2.33 2.15 2.38 -0.452 0.133 0.117 0.0156 
## 3 2.32 2.17 2.43 -1.01 0.122 0.115 0.0705 
## 4 2.80 2.35 2.82 -0.200 0.0582 0.117 0.00124 
## 5 2.93 2.37 2.86 0.650 0.0536 0.116 0.0120 
## 6 2.98 2.38 2.88 0.884 0.0516 0.115 0.0213 
## 7 2.75 2.40 2.92 -1.56 0.0478 0.112 0.0609 
## 8 2.90 2.40 2.92 -0.183 0.0478 0.117 0.000842
## 9 3.12 2.41 2.94 1.57 0.0461 0.112 0.0594 
## 10 2.99 2.42 2.96 0.259 0.0445 0.117 0.00156 
## # … with 21 more rows
```

]

---

# レッスン .circle-big[三]

.tl.ml7.f1[

<ul class="fa-ul">
 <li style="opacity:0.3;">hypothesis testing (frequentist version)</li>
 <li style="opacity:0.3;">linear model</li>
 <li>residual analysis</li>
 <li style="opacity:0.3;">null plots</li>

</ul>
]

---

# Model diagnostic

.aim-box[
`$e_i \sim NID(0, \sigma^2)$` no pattern in **residual plot** & straight line for **QQ-plot** (theoretically)
]

```r
ggplot(res, 
    aes(`log(Girth)`, .std.resid)) + 
  geom_point() + 
  geom_hline(yintercept = 0) + 
  ggtitle("Residual plot")
```

]
.w-50.bl[

```r
ggplot(res, 
    aes(sample = .std.resid)) + 
  geom_qq() +
  geom_qq_line() +
  ggtitle("QQ plot")
```

]

---

]
.w-50.center[
.f2[
And is this line straight in the QQ-plot?
]

]
]

---

<a href="#hypothesis-test" class="button">Hypothesis testing framework</a>

---

## Coin toss example

| | |
|-----|-----|
|**Hypotheses** | `$H_0: p = 0.5$` vs. `$H_1: p > 0.5$` |
|**Assumptions** | Each toss is independent with equal chance of getting a head. |
|**Test statistic** | `$X \sim B(n, p)$` under `$H_0$`. |
|**P-value** .font_small[(or critical value or confidence interval)] | `$P(X \geq x)$` where `$x$` is the observed test statistic. |
|** Conclusion** | Reject null hypothesis when the `$p$`-value is less than some significance level `$\alpha$`. Usually `$\alpha = 0.05$`.|

]

.w-50.bl.pl3[

## Tree example

| | |
|-----|-----|
|**Hypotheses** | `$H_0: e_i \sim NID(0, \sigma^2)$` vs. `$H_1:$` not `$H_0$` |
|**Assumptions** | Assume method of moments estimate for `$\sigma^2$` is good enough. |
|**Test statistic** | Residual and QQ-plot But what is its distribution?? |
|**P-value** | ??? |
|** Conclusion** | ???|

]

---

# レッスン .circle-big[四]

.tl.ml7.f1[

<ul class="fa-ul">
 <li style="opacity:0.3;">hypothesis testing (frequentist version)</li>
 <li style="opacity:0.3;">linear model</li>
 <li style="opacity:0.3;">residual analysis</li>
 <li>null plots</li>

</ul>
]

---

# Samples from the null distribution

]

* The estimate of `$\sigma$`, denoted `$\hat{\sigma}$` is given below:

```r
sigma_hat <- glance(fit)$sigma
sigma_hat
```

```
## [1] 0.1149578
```
{{content}}
--

* Parametric bootstrap `$e_i$`s:

```r
null_generator <- function() {
 rnorm(n = nrow(trees), 
 mean = 0, 
 sd = sigma_hat) 
 }
```

---

# Samples from the null distribution

.flex[
.w-50[
* Under `$H_0$`, we have `$e_i \sim NID(0, \sigma^2)$`.
* The estimate of `$\sigma$`, denoted `$\hat{\sigma}$` is given below:

```r
sigma_hat <- glance(fit)$sigma
sigma_hat
```

```
## [1] 0.1149578
```
* Parametric bootstrap `$e_i$`s:

```r
null_generator <- function() {
 rnorm(n = nrow(trees), 
 mean = 0, 
 sd = sigma_hat) 
 }
```
]
.w-50.f4[

```r
set.seed(2021) # for reproducibility
sample1 <- null_generator()
sample1
```

```
##  [1] -0.014077733  0.063509210  0.040079987  0.041342539  0.103238294 -0.221014401  0.030089561  0.105251514  0.001583192
## [10]  0.198872794 -0.124407911 -0.031363388  0.020921794  0.173418675  0.184446385 -0.211692022  0.186612202  0.015104195
## [19]  0.170266610  0.173967771 -0.108341224 -0.021345944 -0.126582888  0.138882296 -0.186799392  0.012114064 -0.167314594
## [28] -0.040696923 -0.010771552  0.126530466 -0.225757053
```
{{content}}

]

```r
sample2 <- null_generator()
sample2
```

```
##  [1] -0.166452530  0.117192993 -0.163403008 -0.069495694 -0.182032707 -0.147827979 -0.167227403 -0.010009506  0.058023401
## [10]  0.013379792  0.202350333 -0.039673836  0.243710596 -0.003951961 -0.091064303  0.169622012 -0.083408475  0.035910414
## [19]  0.079546685 -0.057512339 -0.259329822  0.005028408 -0.042398523 -0.110385074  0.011928749  0.049120227 -0.019598186
## [28] -0.178085791 -0.173080487  0.001844330 -0.021309077
```

---

.flex.f4[
.w-50[

## Data plot: .monash-blue[Residual plot]

```r
gres <- ggplot(res, 
 aes(`log(Girth)`, .std.resid)) + 
 geom_point() + 
 geom_hline(yintercept = 0) +
 theme_void()

gres
```

<img src="images/unnamed-chunk-12-1.png" width="172.8" style="display: block; margin: auto;" />
 
.center[
Null plot と違いありますか？

]]
.w-50[

## Null plots

```r
gres %+% 
  mutate(res, .std.resid = sample1)
```

```r
gres %+% 
  mutate(res, .std.resid = sample2)
```

<img src="images/unnamed-chunk-14-1.png" width="172.8" style="display: block; margin: auto;" />
]

]

---

.flex.f4[
.w-50[

## Data plot: .monash-blue[QQ-plot]

```r
gres <- ggplot(res, 
 aes(sample = .std.resid)) + 
 geom_qq() + 
 geom_qq_line() +
 coord_equal() +
 theme_void()

gres
```

<img src="images/unnamed-chunk-15-1.png" width="172.8" style="display: block; margin: auto;" />
 
.center[
Null plot と違いありますか？

]]
.w-50[

## Null plots

```r
gres %+% 
  mutate(res, .std.resid = sample1)
```

```r
gres %+% 
  mutate(res, .std.resid = sample2)
```

<img src="images/unnamed-chunk-17-1.png" width="172.8" style="display: block; margin: auto;" />
]

]

---

.tl.ml7.f1[

<ul class="fa-ul">
 <li style="opacity:0.3;">hypothesis testing (frequentist version)</li>
 <li style="opacity:0.3;">linear model</li>
 <li style="opacity:0.3;">residual analysis</li>
 <li style="opacity:0.3;">null plots</li>

</ul>
]

.animated.slideInLeft.f-headline[
**Visual Inference**
]

---

Reading data plots requires calibration

]]

.bottom_abs.f6[
 Buja et al. (2008) “Statistical Inference for Exploratory Data Analysis and Model Diagnostics.” Philosophical Transactions. Series A, Mathematical, Physical, and Engineering Sciences 367 (1906): 4361–83.
 
]

---

.center.f2[
どれか変わってるプロットありますか？
]

---

# Orange Tree Growth

.flex.f4[
.w-50[

```r
ggplot(Orange) +
  aes(age, circumference, color = Tree) + 
  geom_point(size = 3) + 
  geom_smooth(method = "lm") 
```

]
.w-50[
.f3[The fitted model:]

```r
fit <- lm(circumference ~ Tree + age:Tree, 
 data = Orange)
```

]

---

# `nullabor`

```r
library(nullabor)
method <- null_lm(circumference ~ Tree + age:Tree, method = "pboot")
fit <- lm(circumference ~ Tree + age:Tree, data = Orange)
fit_df <- Orange %>% 
 mutate(.resid = residuals(fit))
set.seed(2021)
line_df <- lineup(method, true = fit_df)
```

```
## decrypt("bhMq KJPJ 62 sSQ6P6S2 uD")
```

```r
ggplot(line_df, aes(age, .resid)) + 
  geom_point() + 
  geom_hline(yintercept = 0) +
  facet_wrap(~.sample) + 
  labs(x = "", y = "") +
  theme(axis.text = element_blank(),
        panel.border = element_rect(fill = "transparent"),
        axis.line = element_blank(),
        axis.ticks.length = unit(0, "mm"),
        strip.background = element_rect(fill = "red"),
        strip.text = element_text(color = "white", face = "bold", size = 18))
```

]

---

# Statistical significance of the data plot

.center.f2[
🧑🏻‍⚖️🧑🏼‍⚖️🧑🏽‍⚖️🧑🏾‍⚖️🧑🏿‍⚖️👨🏻‍⚖️👨🏼‍⚖️👨🏽‍⚖️👨🏾‍⚖️👨🏿‍⚖️

多数人の評価をまとめる
]

* Let `$X$` be the number of observers out of `$K$` independent observers picking the data plot from the lineup.

* Assume that `$X \sim B(K, 1/m)$` under `$H_0$`.

* The .monash-blue[**p-value**] of a lineup of size `$m$` evaluated by `$K$` observers is given as `$P(X \geq x)$`.

]]

Unlike conventional hypothesis testing, visual inference p-value depends on

* plot design,
* the individuals' visual skills,

]
.item.pad20[

* the number of `$K$` observers,
* the size `$m$` of the lineup, and

]
.item[
* effect size.
]

]

.bottom_abs.font_small[
Majumder, M., Heiki Hofmann, and Dianne Cook. 2013. “Validation of Visual Statistical Inference, Applied to Linear Models.” Journal of the American Statistical Association 108 (503): 942–56.
]

---

layout: false
background-size: cover
class: title-slide
background-image: url("images/bg-03.png")

This slide is made using the `xaringan` R-package and found at

.w-80.center[
[emitanaka.org/slides/FukuokaR2021](https://emitanaka.org/slides/FukuokaR2021)

]

.bottom_abs.width100[

*Emi Tanaka*

Department of Econometrics and Business Statistics, &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Monash University, Melbourne, Australia

emi.tanaka@monash.edu

<a href="https://twitter.com/statsgen" style="color:black"> @statsgen</a>

]

Thank you!