考古学のためのR統計解析〜Ben Marwick先生の授業を復習する（ANOVA編）〜


この記事は、2020年5月7日に行われた「第2回考古学のためのR統計解析」の復習です。前回の[「カイ2乗検定編」](https://qiita.com/ishiijunpei/items/b27be10622ac8c7a7baa)と[t検定編](https://qiita.com/ishiijunpei/items/d1fb753fb64748ba0cd7)の続編です。

[当日資料「Statistical Inference and Data Exploration for Archaeologists」](https://benmarwick.github.io/stat-inference-and-exploration-for-archaeologists/stat-inference-and-exploration-for-archaeologists.html#1)

使用するデータは、前回と同様、北海道埋蔵文化財センターさんの[厚真町オコッコ1遺跡（２）](https://sitereports.nabunken.go.jp/70356)（北埋調報356,p234-243）の石器一覧表データを使用しました。`tabulizer`パッケージを使用して、pdfの表データを読み込んでいます。

```{r}
library(tidyverse)
library(infer)
library(knitr)
```

```{r}
# データ読込
tbs <-
  read_csv("myData/table2.csv")
```

長さと幅を実数型に変換しつ、つまみ付きナイフを抽出します。

```{r}
tbs4 <-
  tbs %>%
  mutate(
    長さ = as.numeric(長さ) ,
    幅 = as.numeric(幅) ,
    ) %>%
  filter(遺物名 %in% "つまみ付きナイフ")
```

前回のt検定は、2つのカテゴリー間における、連続量の差の有無を確認するものでした。ANOVA（分散分析）は3つ以上のカテゴリー間において差の有無とどのカテゴリー間に差があるのかを確認します。分析対象はまたしても「つまみ付きナイフ」です。石材ごとに長さに差があるかどうかを見ていきます。

```{r}
# boxplotで可視化
tbs4 %>%
  ggplot(aes(x = 石材 , y = 長さ))  +
  geom_boxplot()

```

![boxplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/6a18c0d9-4fea-68a5-4bf2-749223dd5a19.png)

Sh(頁岩)が他の石材より若干長いように見えることは、t検定でも確認済みです。Obs(黒曜石)は短いように思えます。

おなじみの手順です。F値を算出してf_statに付値します。

```{r}
# compute F statistic
f_stat <-
  tbs4 %>%  
    specify(長さ ~ 石材) %>%  # 長さと石材の列をtibble形式で取り出します。
    calculate(stat = "F")  # F値を算出します。
```
リサンプリングによって、石材と長さの組み合わせを1000パターン抽出します。

```{r}
# 1000回のリサンプリングでF値を算出
null_distribution <-  
  tbs4 %>%
    specify(長さ ~ 石材) %>% # 長さと石材の2列のtbble形式データを作成
    hypothesize(null = "independence") %>% # 帰無仮説は「石材と長さは独立である」
    generate(reps = 1000, type = "permute") %>% # 1000回リサンプリング
    calculate(stat = "F") # F値を算出
```

```{r}
# 可視化
null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = f_stat , direction = "right")
```
![histogram.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/eaa81485-6028-ebd7-8e62-e6f00bc126eb.png)

これは、帰無仮説が棄却されない予感がします。

```{r}
aov(長さ ~ 石材,  
    data = tbs4) %>%  
  broom::tidy() %>%　# リスト形式のデータをtibble形式に変換
  kable()
```
p値は0.14ですから、帰無仮説は棄却されません。

|term      | df|     sumsq|   meansq| statistic|   p.value|
|:---------|--:|---------:|--------:|---------:|---------:|
|石材      |  3|  17.67105| 5.890349|  1.820653| 0.1483963|
|Residuals | 98| 317.05885| 3.235294|        NA|        NA|


この様子を石材の組み合わせごとに可視化します。

```{r}
# 可視化
aov(長さ ~ 石材,  
    data = tbs4) %>%  # 分散分析を行うaovを実行
  TukeyHSD() %>%  # チューキー・クレーマー検定を実行
  broom::tidy() %>%  #TukeyHSD独自のクラスをtibble形式に変換。リスト形式をテーブル形式に直している。
  ggplot(aes(x = fct_reorder(comparison,  # 石材の組み合わせをestimate(平均値の差)の昇順に整列
                             estimate),  
             y = estimate)) +  # 石材同士の平均値の差
  geom_pointrange(aes(ymin = conf.low,  # 信頼区間の下限値
                      ymax = conf.high,  # 信頼区間の上限値
                      color = adj.p.value < 0.05)) +  # p値が0.05以下で色を分ける
  geom_hline(yintercept = 0,  # y切片0の直線
             linetype = "dashed") +  # 点線
  coord_flip() +
  labs(x = NULL,  
       colour = "padj < 0.05") 
```
![ANOVAplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e6eaa4f7-d780-4ff5-5a49-d6807d957a18.png)

赤丸は石材の組み合わせにおける平均値の差、エラーバーの範囲にestimate=0の点線が含まれている組み合わせは、信頼区間内に0（差がない）が含まれています。
グラフを解釈するために、ggplotへの入力を見ておきます。

```{r}
aov(長さ ~ 石材,  
    data = tbs4) %>%  # 分散分析を行うaovを実行
  TukeyHSD() %>%  # チューキー・クレーマー検定を実行
  broom::tidy() %>% 
  kable()
```

|term |comparison |   estimate|   conf.low| conf.high| adj.p.value|
|:----|:----------|----------:|----------:|---------:|-----------:|
|石材 |Obs-Ag     | -0.4714286| -4.2407575| 3.2979004|   0.9878671|
|石材 |Sh-Ag      |  0.6917808| -2.6776846| 4.0612463|   0.9498943|
|石材 |Si-Ag      | -0.1450000| -3.6314878| 3.3414878|   0.9995373|
|石材 |Sh-Obs     |  1.1632094| -0.6969123| 3.0233310|   0.3642824|
|石材 |Si-Obs     |  0.3264286| -1.7381179| 2.3909751|   0.9760774|
|石材 |Si-Sh      | -0.8367808| -2.0232920| 0.3497304|   0.2595695|

- comparison：石材の組み合わせ
- estimate：石材同士の平均値の差
- conf.low：信頼区間の下限
- conf.high：信頼区間の上限
- adj.p.value：p値

estimateが0から離れるほど2つの石材の平均値に差があります。
信頼区間の上限と下限の間に0がなければ、平均値に差がない状態である0は、当該信頼区間の外にあるということになります。
上記のグラフでは、いずれの組み合わせにおいても95%信頼区間に0を含んでいることから、石材によって石器の長さに差があるという帰無仮説は、5%の有意確率で棄却することはできません。

以上、ANOVAでした。難しかった・・・

（次回は「主成分分析編」です。）