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

この記事は、2020年5月7日に行われた「第2回考古学のためのR統計解析」の復習です。前回の「カイ2乗検定編」の続編です。

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

当日資料「Statistical Inference and Data Exploration for Archaeologists」

# パッケージ読み込み
library(tidyverse)
library(infer)
library(knitr)

データ準備

1. 長さと幅を実数型に変換して「つまみ付きナイフ」だけを選択します。
2. 石材をShとそれ以外（Other）に分けます。

tbs2 <-
  tbs %>%
  mutate(
    長さ = as.numeric(長さ) , # 長さを実数に変換
    幅 = as.numeric(幅) , # 幅を実数に変換
    ) %>%
  filter(遺物名 %in% "つまみ付きナイフ") %>% # つまみ付きナイフだけを選択
  mutate( 石材分類 = # mutate()を利用して石材を再分類（頁岩とそれ以外にする）
    case_when(
      石材 %in% "Sh" ~ "Sh" , # 石材　sh（頁岩）はshにする（そのまま）
      石材 %in% c("Ag" , "Obs" , "Si") ~ "Other" # Ag Obs Si（メノウ、黒曜石、珪岩？） をOtherにする。
    )
  )

# クロス集計表
tbs2 %>%
  count(遺物名, 石材分類) %>%
  spread(key = 石材分類 ,value = n) %>%
  kable()

作成したデータはこのような形です。Sh（頁岩）が圧倒的に多いです。

遺物名	Other	Sh
つまみ付きナイフ	29	73

基本的な手順はカイ2乗検定と同じす。標本から1000回リサンプリングした標本群に対して、t値を算出します。

# 1000回のリサンプリングでt値を算出
null_distribution <-  
  tbs2 %>%
  specify(長さ ~ 石材分類) %>% # 長さと石材分類の2列のtbble形式データを作成
  hypothesize(null = "independence") %>% # 帰無仮説は「石材分類と長さは独立である」
  generate(reps = 1000, type = "permute") %>% # 1000回リサンプリング
  calculate(stat = "t", order = c("Sh", "Other")) # t統計量を計算

続いて、データ本体のt値を算出し、t_statオブジェクトに格納します。

# 標本のt値を算出
t_stat <- 
  tbs2 %>%
  specify(長さ ~ 石材分類) %>% 
  calculate(stat = "t" ,order = c("Sh", "Other")) %>%
  dplyr::pull() # caluculate()の結果の要素（chitestの統計量）にアクセス

infer::visualize()でt値のヒストグラムを表示します。

# 可視化する
null_distribution %>% # 1000回サンプリングした標本のt値
  visualize() + # ヒストグラム描画
  shade_p_value(obs_stat = t_stat, # 標本のt値
                direction = "right") # 片側検定

見た目には、帰無仮説が棄却されるかどうか、なかなか微妙なところにいます。
帰無仮説は棄却されるでしょうか？

密度図で石材分類ごとの長さの分布を比較します。

# 密度図で可視化する
tbs2 %>%
  ggplot(aes(x = 長さ , fill = 石材分類)) +
  geom_density(alpha = 0.4 , colour = "white") +
  scale_fill_viridis_d(option = "cividis")

密度図からは、石材によって石器の長さが明らかに異なっている、という印象は受けません。

いよいよp値を算出します。

# p値の算出
tbs2 %>%
  t_test(長さ ~ 石材分類 , order = c("Sh","Other")) %>%
  kable()

statistic	t_df	p_value	alternative	lower_ci	upper_ci
-2.589326	66.56431	0.0118003	two.sided	-1.603729	-0.2074192

「棄却すべき有意確率は最初に決めておきなさい」という教えの理由がよくわかる結果です。5%なら棄却できるけど、1％では棄却できません。

（次回は「ANOVA編」です。）