이 기사는 2020년 5월 7일에 열린 '제2회 고고학을 위한 R 통계 분석'의 복습입니다. 마지막 「카이 제곱 검정편」의 속편입니다.

사용하는 데이터는, 전회와 같이, 홋카이도 매장 문화재 센터씨의 아츠 마마치 오코코 1 유적 (2) (북매조보 356, p234-243)의 석기 일람표 데이터를 사용했습니다. tabulizer 패키지를 사용하여 pdf 테이블 데이터를 읽는 중입니다.

당일 자료 「Statistical Inference and Data Exploration for Archaeologists」

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

데이터 준비

길이와 너비를 실수형으로 변환하고 노브 나이프만 선택합니다.

석재를 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(혈암)가 압도적으로 많습니다.


유물명
기타
Sh


노브가 달린 칼
29
73


기본 절차는 카이 제곱 검정과 동일합니다. 표본에서 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편」입니다.)