この記事は、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検定編です（先は長い）。