この記事は、2020年5月7日に行われた「第2回考古学のためのR統計解析」の復習です。前回の「カイ2乗検定編」とt検定編の続編です。
当日資料「Statistical Inference and Data Exploration for Archaeologists」
使用するデータは、前回と同様、北海道埋蔵文化財センターさんの厚真町オコッコ1遺跡（２）（北埋調報356,p234-243）の石器一覧表データを使用しました。tabulizerパッケージを使用して、pdfの表データを読み込んでいます。
library(tidyverse)
library(infer)
library(knitr)
# データ読込
tbs <-
read_csv("myData/table2.csv")
長さと幅を実数型に変換しつ、つまみ付きナイフを抽出します。
tbs4 <-
tbs %>%
mutate(
長さ = as.numeric(長さ) ,
幅 = as.numeric(幅) ,
) %>%
filter(遺物名 %in% "つまみ付きナイフ")
前回のt検定は、2つのカテゴリー間における、連続量の差の有無を確認するものでした。ANOVA（分散分析）は3つ以上のカテゴリー間において差の有無とどのカテゴリー間に差があるのかを確認します。分析対象はまたしても「つまみ付きナイフ」です。石材ごとに長さに差があるかどうかを見ていきます。
# boxplotで可視化
tbs4 %>%
ggplot(aes(x = 石材 , y = 長さ)) +
geom_boxplot()
Sh(頁岩)が他の石材より若干長いように見えることは、t検定でも確認済みです。Obs(黒曜石)は短いように思えます。
おなじみの手順です。F値を算出してf_statに付値します。
# compute F statistic
f_stat <-
tbs4 %>%
specify(長さ ~ 石材) %>% # 長さと石材の列をtibble形式で取り出します。
calculate(stat = "F") # F値を算出します。
リサンプリングによって、石材と長さの組み合わせを1000パターン抽出します。
# 1000回のリサンプリングでF値を算出
null_distribution <-
tbs4 %>%
specify(長さ ~ 石材) %>% # 長さと石材の2列のtbble形式データを作成
hypothesize(null = "independence") %>% # 帰無仮説は「石材と長さは独立である」
generate(reps = 1000, type = "permute") %>% # 1000回リサンプリング
calculate(stat = "F") # F値を算出
# 可視化
null_distribution %>%
visualize() +
shade_p_value(obs_stat = f_stat , direction = "right")
これは、帰無仮説が棄却されない予感がします。
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 |
この様子を石材の組み合わせごとに可視化します。
# 可視化
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")
赤丸は石材の組み合わせにおける平均値の差、エラーバーの範囲にestimate=0の点線が含まれている組み合わせは、信頼区間内に0（差がない）が含まれています。
グラフを解釈するために、ggplotへの入力を見ておきます。
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でした。難しかった・・・
（次回は「主成分分析編」です。）