고고학을 위한 R 통계 분석 ~ Ben Marwick 선생님의 수업을 복습한다(ANOVA편)~

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

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

사용하는 데이터는, 전회와 같이, 홋카이도 매장 문화재 센터씨의 아츠 마마치 오코코 1 유적 (2) (북매조보 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-검정은 두 범주 간의 연속량 차이의 유무를 확인하는 것이었다. 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에서 멀어질수록 두 석재의 평균값에 차이가 있습니다.
신뢰 구간의 상한과 하한 사이에 0이 없으면 평균값에 차이가 없는 상태인 0은 해당 신뢰 구간 밖에 있다고 합니다.
위의 그래프에서는 어떠한 조합에 있어서도 95% 신뢰 구간에 0을 포함하고 있기 때문에, 석재에 의해 석기의 길이에 차이가 있다고 하는 귀무 가설은, 5%의 유의 확률로 기각할 수 없습니다 .

이상, ANOVA였습니다. 어려웠다・・・

(다음 번은 「주성분 분석편」입니다.)