この記事は、2020年5月7日に行われた「第2回考古学のためのR統計解析」の復習です。
当日は、ものすごスピードで盛りだくさんの内容が詰め込まれていたので、ゆっくりと勉強してみました。
当日資料「Statistical Inference and Data Exploration for Archaeologists」
パッケージ読み込み
library(tidyverse)
報告書データ読み込み
セミナーでは、Marwick先生がデータを用意されていましたが、せっかくなので日本の北海道のデータを使用します。
北海道埋蔵文化財センターさんの厚真町オコッコ1遺跡（２）（北埋調報356,p234-243）の石器一覧表データを使用しました。
- 全国遺跡総覧から報告書データをダウンロード
- pdftkコマンドで一覧表部分を抽出（Rコードではありません）
#247-256ページを抽出
pdftk myData/okoko.pdf cat 247-256 output myData/table.pdf
一覧表PDFはこんな感じです。
tabulizerパッケージでpdfのtableデータをtibble形式に変換し、csvに書き出します。
library(tabulizer)
# pdfデータをtibble形式に変換
tabulizer::extract_tables("myData/table.pdf") %>%
purrr::map_dfr(as.data.frame) %>%
as_tibble()%>%
write_csv("myData/table.csv")
Rだけでデータクリーニングができなかったので、リブレオフィスでcsvを開いて目視で編集します。
この時点で、残念ながら再現性が破綻してしまいました。
# 手動でデータ編集して再読込（再現性破綻）
tbs <-
read_csv("myData/table2.csv")
# 読み込んだデータ
> tbs
# A tibble: 762 x 12
遺構 発掘区 小グリッド 層位 取り上げ日 遺物名 細分記号 石材 長さ 厚さ 幅
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <chr>
1 北盛土… K A MII上 2015/8/21 石槍 Ib Sh 5.8 0.7 2.2
2 南盛土… T C MII上 2015/9/18 石槍 Ib Obs 8 1.1 2.9
3 NA M NA V 2015/7/6 石槍 Ib Sh 7.8 0.9 3.1
4 北盛土… J D MII中 2015/9/10 石槍 Ib Sh 8 0.8 3.1
5 南盛土… S NA MI 2015/8/7 石槍 Ib Obs 8.3 1.3 3.3
6 南盛土… R NA MI 2015/8/28 石槍 Ib Sh 6.3 1 3.2
7 NA M NA V上 2015/8/4 石槍 Ib Sh 8.4 1 2.8
8 NA S NA NA 2015/8/17 石槍 Ib Sh 10.60 0.8 3.8
9 NA T NA V 2015/7/8 石槍 Ia Sh 6.7 1 2.7
10 NA O NA V中 2015/10/1 石槍 Ib Obs 5.4 0.7 2.5
# … with 752 more rows, and 1 more variable: 重さ <chr>
カイ2乗検定
準備として、剝片石器を3種類選択します。オーソドックスな3器種を選びました。
# 剝片石器をいくつか選択
tbs3 <-
tbs2 %>%
filter(遺物名 %in% c("石槍", "石鏃" ,"スクレイパー"))
library(infer)は仮説検定のためのモダンなパッケージです。ggplot2と連携して、仮説検定を視覚的に把握することができます。以下で行っているのは次のような手順です。
-
specify()で分析対象の「遺物名」と「石材」のみのtableを生成しつつ、factor型に変換します。 -
hypothesize()で、帰無仮説を指定します。ここでは「遺物名と石材は独立である」が、帰無仮説です。 -
generate()で、標本からのリサンプリングを行います。今回は1000回のリサンプリングを指定しました。 -
caluculate()のオプションstat = に"Chisq"を指定します。
library(infer) # 統計的検定のためのモダンなパッケージ
null_distribution <-
tbs3 %>%
specify(遺物名 ~ 石材) %>% # 遺物名と石材をfactor型に変換
hypothesize(null = "independence") %>% # 独立性の検定を指定
generate(reps = 1000, type = "permute") %>% # reps= でリサンプリング回数を指定 type = "bootstrap", "permute"," or "simulate" が選択可能
calculate(stat = "Chisq") # カイ2乗統計量を算出
標本本体のカイ2乗検定の統計量を算出します。
chi_square_stat <-
tbs3 %>%
specify(遺物名 ~ 石材) %>%
calculate(stat = "Chisq") %>%
dplyr::pull() # caluculate()の結果の要素（chitestの統計量）にアクセス
先ほど1000回リサンプリングした標本の分布にカイ2乗検定の統計量（chi_square_stat=2.88ぐらい）を重ねます。
標本分布の真ん中ほどに赤いバーで示されているのが、データから直接算出されたカイ2乗の統計量です。
このグラフを見る限りでは、「遺物名と石材は独立である」という帰無仮説を棄却することは難しそうです。
null_distribution %>% # 1000回リサンプリングしたカイ2乗の統計量
visualize() + # ヒストグラム描画
shade_p_value(obs_stat = chi_square_stat, # データの統計量
direction = "right") # 片側（上側）検定を指定
p値を算出してみます。
library(knitr)
tbs3 %>%
chisq_test(遺物名 ~ 石材) %>% # 遺物名と石材の独立性の検定
kable() # 表に書き出し
| statistic | chisq_df | p_value |
|---|---|---|
| 2.889206 | 3 | 0.4090247 |
p_value = 0.4090 ですから、やはり帰無仮説を棄却することはできなさそうです。
つまり、「石器の器種によって石材の偏りがあると判断することはできない」という結論になります。
次回はt検定編です（先は長い）。