ggplot2のグラフ表現を改良する

奈良文化財研究所で9月10日〜11日に開催された[「考古学・文化財のためのデータサイエンス研究集会」](https://sitereports.nabunken.go.jp/ja/event/332)で発表した資料[考古学のためのデータビジュアライゼーション](https://github.com/IshiiJunpei/2019datasience)（ggplot2を使用したデータ可視化）について、[Ben Marwick](https://twitter.com/benmarwick)先生から指導をいただきましたので、備忘録として公開します。


# alphaはaes()の外に出す
下記のコードでは`alpha＝0.7`に設定していますが、alphaの指定を`aes()`の中で行っているので、無意味な凡例がついてしまっています。
あと、`pirnt(p)`は明示しなくても良いのでは？という提案もいただきました。

```{r}
p<-pot%>%
	ggplot(aes(x=口径,fill=分類,alpha=0.7))+
		geom_density()+
		scale_fill_ptol()+
		theme_minimal()
print(p)
```
![density_befor.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/ccbe71cb-5f3d-cff3-3efe-f1b1126c2357.png)


修正したのがこれ。`alpha`を`aes()`の中に入れるのではなく、`geom_density()`の引数として指定しました。コードもスッキリしています。

```{r}
pot%>%
	ggplot(aes(x=口径,fill=分類))+
		geom_density(alpha=0.7)+
		scale_fill_ptol()+
		theme_minimal()
```
![density_after.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/81ae8f12-64f1-a592-a0ed-e5825561a5ef.png)

# 箱ひげ図を改良する

次の指摘はこちらの箱ひげ図。Marwick先生のご指摘は2点。

- ggforceパッケージの`geom_sina()`を使って箱ひげ図とドットを重ねて表示する。
- カラーパレットは [viridis colour palettes](https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html)を使ったほうが良い。

ggforceパッケージは初めて知りました。viridis colourは連続量のイメージがありますが、離散的なデータにも対応しているようです。

```{r}
pot%>%
	ggplot(aes(x=分類,y=口径,fill=分類))+
		geom_boxplot()+
		scale_fill_ptol()+
		coord_flip()+
		theme_minimal()
```

![boxplot_bifore.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/f30feae5-6b2d-075e-bf03-a89d9f13442b.png)

ということで改良してみました。
viridisカラーパレットを使用する理由として、単に見やすいということだけではなく、ユニバーサルなカラーパレットであるということが挙げられます。brewerやptolも優れたカラーパレットですが、viridisカラーが現時点では一日の長があるようです。


```{r}
pot%>%
	ggplot(aes(x=分類,y=口径,fill=分類))+	
		geom_boxplot(alpha = 0.2)+		#不透明度を0.2
    	geom_sina(aes(colour = 分類),		#geom_sina()関数でaes()の引数にcolour=分類を指定
              alpha = 0.4, 
              size = 3) +
		scale_fill_viridis_d() +		#viridis_d()は離散量、連続量ならviridis_c()を指定する
  		scale_colour_viridis_d() +
		coord_flip()+
		theme_minimal()
```
![boxplot_after.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/53e297f9-92b3-c2eb-55f7-0aaa0204d061.png)

# 統計データは小数点2〜3桁で表示すれば十分

分散分析の結果を表示した表に対してのコメントです。こうした統計表の表示には小数点2〜3桁で十分とのこと。

```{r}
# aov関数の結果をTukeyHSD関数に渡す
tkh<-aov(口径~分類,data=pot)%>%TukeyHSD()
tkh$分類%>%kable(format="markdown")
```

<table>
<thead>
<tr class="header">
<th style="text-align: left;"></th>
<th style="text-align: right;">diff</th>
<th style="text-align: right;">lwr</th>
<th style="text-align: right;">upr</th>
<th style="text-align: right;">p adj</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td style="text-align: left;">B型-A型</td>
<td style="text-align: right;">-6.58</td>
<td style="text-align: right;">-8.1885528</td>
<td style="text-align: right;">-4.971447</td>
<td style="text-align: right;">0.0000000</td>
</tr>
<tr class="even">
<td style="text-align: left;">C型-A型</td>
<td style="text-align: right;">-4.54</td>
<td style="text-align: right;">-6.1485528</td>
<td style="text-align: right;">-2.931447</td>
<td style="text-align: right;">0.0000000</td>
</tr>
<tr class="odd">
<td style="text-align: left;">C型-B型</td>
<td style="text-align: right;">2.04</td>
<td style="text-align: right;">0.4314472</td>
<td style="text-align: right;">3.648553</td>
<td style="text-align: right;">0.0087802</td>
</tr>
</tbody>
</table>

改良したのがこちら。小数点3桁でまるめられて良い感じ。ただ、もともとあった変数の組み合わせ（"B型-A型"の部分）が消えてしまっているのが少し残念です。

```{r}
tkh <- 
	aov(口径 ~ 分類, data = pot) %>% 
	TukeyHSD() %>% 
	.$分類 %>%	#TukeyHSD関数の結果から$分類を選択 
	as_tibble() %>% 	#tibble_df形式に変換
	mutate_if(is.numeric, round,3)		#mutate_if()でnumericクラスのカラムにround関数を適用する。
tkh%>%kable(format="markdown")
```

|  diff|    lwr|    upr| p adj|
|-----:|------:|------:|-----:|
| -6.58| -8.189| -4.971| 0.000|
| -4.54| -6.149| -2.931| 0.000|
|  2.04|  0.431|  3.649| 0.009|

# 集計表タイプのデータはgeom_barのかわりにgeom_colを使う

`geom_bar(stat="identity")`を指定していますが、`geom_col(position="fill")`とするべき、との提案です。結果はどちらも同じですが、`geom_col`は`stat="identity"`がデフォルトとなっていますので、集計表タイプのデータ描画に利用しやすくなっています。

befor

```{r before}
toj%>%
	ggplot(aes(x=遺跡名,y=点数,fill=器種))+
		geom_bar(stat="identity",position="fill")+
		coord_flip()+
		scale_fill_ptol()+
		theme_minimal()
```
after

```{r}
toj%>%
	ggplot(aes(x=遺跡名,y=点数,fill=器種))+
		geom_bar(stat="identity",position="fill")+
		coord_flip()+
		scale_fill_ptol()+
		theme_minimal()
```

![bar_before.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/5772bdc5-e49e-9a48-31e8-f636e9f1358a.png)

# tidyverseだけでデータ整形

細別器種を大別器種に再分類していますが、コードが冗長です。Marwick先生は次のように指摘します。

- `grepl()`を多用しているが、もっとシンプルに記述できる。
- `tidyverse`で全てを処理したいなら`glepl()`ではなく`stringr::str_detect()`を使っても良い。


```{r}
toj2<-toj%>%
	mutate(大別器種 = case_when(
		grepl("碗",器種)|grepl("皿",器種)|grepl("その他の食膳具",器種)	 == TRUE ~ "食膳具",
		grepl("壺･甕･瓶",器種) == TRUE ~ "貯蔵具",
		grepl("灯明皿・油注",器種)|grepl("その他",器種)|grepl("すり鉢",器種)|
			grepl("鍋",器種)|grepl("土瓶",器種) == TRUE ~ "その他",
		))
```

次のように修正します。シンプル。

```{r}
toj2<-toj%>%
	mutate(
		大別器種 = case_when(
			str_detect(器種,"碗 | 皿 | その他食膳具") ~ "食膳具",
			str_detect(器種,"壺･甕･瓶") ~ "貯蔵具",
			str_detect(器種,"灯明皿・油注 | その他 | すり鉢 | 鍋 | 土瓶") ~ "その他",
		)
	)
```
大別器種というカラムが作られます。

|遺跡名         |器種         | 点数|大別器種 |
|:--------------|:------------|----:|:--------|
|弁天貝塚       |碗           |  134|食膳具   |
|弁天貝塚       |皿           |   84|食膳具   |
|弁天貝塚       |その他食膳具 |   34|食膳具   |
|弁天貝塚       |土瓶         |    6|その他   |
|弁天貝塚       |鍋           |    0|その他   |

# 散布図行列に有意性の評価を入れる

`GGally::ggpairs()`は見栄えのよい表現ですが、有意性の評価を図中に入れるため`PerformanceAnalytics::chart.Correlation()`を使用してはどうか、という提案をいただきました。

`GGally::ggpairs()`による描画
![pairs_before.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/ca29c3c3-bccf-0dd5-86c8-7a8d8bc0547c.png)

`PerformanceAnalytics::chart.Correlation()`による描画です。文字サイズの調整のやり方などがうまくできていませんが、有意性の評価が星印で示されています。

```{r}
library(PerformanceAnalytics)
iron %>%
	select(全長, 刀身長, 茎長, 刀身先幅, 刀身元幅, 刀身元厚, 茎先幅) %>%
		chart.Correlation(histogram = TRUE, pch = 19)
```
![pairs_after.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/d8c04780-df70-a3b2-e89f-79cef9337b4e.png)

# 図中に数式を描画する

研修では以下のように`lm()`を適用して必要なパラメーターを算出して手動で回帰式を表示していました。

```{r}
icoe<-lm(刀身長 ~ 刀身元幅,data=iron)%>%summary()
icoe$coefficients%>%kable()
```
|            |  Estimate| Std. Error|   t value| Pr(>&#124;t&#124;)|
|:-----------|---------:|----------:|---------:|------------------:|
|(Intercept) | -6.280888|  1.9780720| -3.175257|          0.0022892|
|刀身元幅    | 10.723991|  0.7889999| 13.591878|          0.0000000|

y=10.72399X-6.280888

`ggpmisc::stat_poly_eq()`と`ggpmisc::stat_fit_glance()`を利用して回帰式、相関係数、p値をプロットに記述する手法を教えていただきました。
複雑な記述になるので、引数やパラメーターの意味が理解できていませんが、再現性を保つためにはこの手法が欠かせないところです。

```{r}
library(ggpmisc)
iron %>%
	ggplot(aes(x=刀身元幅,y=刀身長))+
		geom_point()+
		geom_smooth(method="lm")+
		theme_minimal() +
		stat_poly_eq(formula = y ~ x,
			eq.with.lhs = "italic(hat(y))~`=`~",
			aes(label = paste(stat(eq.label), 
				stat(rr.label), 
				sep = "~~~")
			), parse = TRUE
		) +
		stat_fit_glance(label.y = 0.9,
			method = "lm",
			method.args = list(formula = y ~ x),
			aes(label = sprintf(
				'~~italic(P)~"="~~%.25f',
				stat(p.value)
				)
			),parse = TRUE
		)
```
![lm_after.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/23875aff-0a84-61d6-8dd6-f540feca1e66.png)