この記事は、2020年5月7日に行われた「第2回考古学のためのR統計解析」の復習です。前回の「カイ2乗検定編」とt検定編と「ANOVA編」の続編です。
当日資料「Statistical Inference and Data Exploration for Archaeologists」
使用するデータは、前回と同様、北海道埋蔵文化財センターさんの厚真町オコッコ1遺跡（２）（北埋調報356,p234-243）の石器一覧表データを使用しました。tabulizerパッケージを使用して、pdfの表データを読み込んでいます。
library(tidyverse)
library(knitr)
厚真町オコッコ1遺跡の石器データを読み込みます。
# データ読込
tbs <-
read_csv("myData/table2.csv")
必要な石器のカテゴリと変数だけを抽出して、実数型に変換します。
tbs5 <-
tbs %>%
filter(遺物名 %in% "つまみ付きナイフ") %>% # つまみ付きナイフを抽出
select(長さ,厚さ,幅,重さ) %>% # 数値属性だけを抽出
mutate_all(funs(as.numeric)) # 全ての変数を実数値に変換
主成分分析を行うモダンなパッケージFactoMineRを読み込み、主成分分析を実行します。
# ライブラリ読み込み
library(FactoMineR)
# 主成分分析を実行
res.pca <-
PCA(tbs5,graph = FALSE)
ここから、主成分分析の結果を使った可視化を行います。
まずはスクリープロットを描画します。
# 多変量解析のヴィジュアライズに特化したfactoextraパッケージ
library(factoextra)
# 各主成分の寄与率を描画します。
fviz_screeplot(res.pca)
第一主成分が60%近くを占めています。第2主成分まで合わせると80%を超えるでしょうか。主成分分析の結果を詳しく調べてみます。
summary(res.pca)
Call:
PCA(X = tbs5, graph = FALSE)
Eigenvalues
Dim.1 Dim.2 Dim.3 Dim.4
Variance 2.296 1.083 0.531 0.089
% of var. 57.410 27.085 13.286 2.218
Cumulative % of var. 57.410 84.496 97.782 100.000
Individuals (the 10 first)
Dist Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3 ctr cos2
1 | 2.107 | 1.813 1.403 0.740 | -0.982 0.873 0.217 | 0.042 0.003 0.000 |
2 | 2.730 | 2.276 2.211 0.695 | -0.301 0.082 0.012 | 1.472 3.996 0.291 |
3 | 2.868 | 2.403 2.465 0.702 | -1.045 0.989 0.133 | 1.165 2.506 0.165 |
4 | 2.679 | 2.005 1.716 0.560 | -1.003 0.911 0.140 | 1.443 3.840 0.290 |
5 | 2.642 | 1.618 1.118 0.375 | -1.340 1.624 0.257 | 1.586 4.642 0.361 |
6 | 1.100 | -0.899 0.345 0.668 | -0.613 0.340 0.310 | 0.160 0.047 0.021 |
7 | 1.525 | -0.736 0.231 0.233 | -0.756 0.517 0.245 | -1.098 2.222 0.518 |
8 | 1.325 | 1.163 0.578 0.770 | 0.499 0.225 0.142 | 0.160 0.047 0.014 |
9 | 3.360 | 3.313 4.685 0.972 | -0.164 0.024 0.002 | 0.529 0.517 0.025 |
10 | 2.730 | -2.169 2.009 0.631 | -1.409 1.796 0.266 | -0.846 1.320 0.096 |
Variables
Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3 ctr cos2
長さ | 0.785 26.846 0.617 | -0.436 17.558 0.190 | 0.416 32.513 0.173 |
厚さ | 0.822 29.459 0.677 | 0.019 0.034 0.000 | -0.559 58.855 0.313 |
幅 | 0.261 2.978 0.068 | 0.941 81.769 0.886 | 0.197 7.271 0.039 |
重さ | 0.967 40.717 0.935 | 0.083 0.639 0.007 | 0.085 1.360 0.007 |
主成分分析の概要はEigenvaluesという要素に格納されているようなので、res.pca$eigを表示します。
res.pca$eig %>%
kable()
| eigenvalue | percentage of variance | cumulative percentage of variance | |
|---|---|---|---|
| comp 1 | 2.2964165 | 57.410413 | 57.41041 |
| comp 2 | 1.0834152 | 27.085381 | 84.49579 |
| comp 3 | 0.5314524 | 13.286311 | 97.78210 |
| comp 4 | 0.0887158 | 2.217896 | 100.00000 |
eigenvaluesは主成分の分散、 percentage of variancev は寄与率、cumulative percentage of variance が累積寄与率を示しています。スクリープロットを作成するfviz_screeplot()は、自動的にpercentage of varianceをy値に出力します。
続いて、主成分1に対する各変数の寄与率を出図します。
fviz_contrib(res.pca,
choice = "var", # 変数ごとの寄与率(ctr)
axes = 1, # 主成分1を指定
top = 10) # 表示する変数の数を指定
第１主成分は重さ、厚さ、長さが高い寄与率を占めていることから、第１主成分は「大きさ」と要約できそうです。
y軸に指定されている"var"でres.pcaオブジェクトの要素であるres.pca$varを引数に指定しています。実際に確認します。
res.pca$var$contrib %>%
kable()
| Dim.1 | Dim.2 | Dim.3 | Dim.4 | |
|---|---|---|---|---|
| 長さ | 26.846489 | 17.5576298 | 32.513332 | 23.082550 |
| 厚さ | 29.459089 | 0.0344304 | 58.855124 | 11.651356 |
| 幅 | 2.977506 | 81.7690726 | 7.271427 | 7.981994 |
| 重さ | 40.716916 | 0.6388672 | 1.360117 | 57.284100 |
各主成分における変数ごとの寄与率がres.pca$var$contribに格納されていることがわかります。Dim.1を棒グラフに出力したものがfviz_contrib()の出力です。
主成分得点の散布図を出力します。
fviz_pca_biplot(res.pca)
右側にプロットされる個体は「重さ」、「厚さ」、「長さ」が大きな値をとる、すなわち「大きい」個体です。第2主成分は「幅」と相関を持ちますので、上にプロットされる個体は「幅広」の個体です。
91〜100ぐらいの個体が少し離れて左上に分布しているのが気になりますね。
主成分1と3を表示する場合は、axes =で指定します。
fviz_pca_biplot(res.pca ,
axes = c(1,3)
)
石材ごとに主成分得点をビジュアル化する
せっかくなので、石材ごとの属性がどのように分布するかを見てみたいです。そのためには、主成分得点が格納されている場所を突き止めます。
res.pca
name description
1 "$eig" "eigenvalues"
2 "$var" "results for the variables"
3 "$var$coord" "coord. for the variables"
4 "$var$cor" "correlations variables - dimensions"
5 "$var$cos2" "cos2 for the variables"
6 "$var$contrib" "contributions of the variables"
7 "$ind" "results for the individuals"
8 "$ind$coord" "coord. for the individuals"
9 "$ind$cos2" "cos2 for the individuals"
10 "$ind$contrib" "contributions of the individuals"
11 "$call" "summary statistics"
12 "$call$centre" "mean of the variables"
13 "$call$ecart.type" "standard error of the variables"
14 "$call$row.w" "weights for the individuals"
15 "$call$col.w" "weights for the variables"
res.pcaに格納された主成分分析結果のうち、主成分得点が格納されているのはres.pca$ind$coordっぽいです。
res.pcaから主成分得点を取り出して、元の石器データと結合します。
# 主成分得点を取り出す
dim1 <- res.pca$ind$coord[,1] # 主成分1の主成分得点
dim2 <- res.pca$ind$coord[,2] # 主成分2の主成分得点
# 石材変数を取り出す
material <-
tbs %>%
filter(遺物名 %in% "つまみ付きナイフ") %>% # つまみ付きナイフを抽出
dplyr::select(石材) # 石材のカラムを抽出
# 石器のデータに付値
material2 <-
material %>%
bind_cols(dim1) %>% # dim1を結合
rename(Dim1 = `NA`) %>% # カラム名をDim1に変更（デフォルトはNA）
bind_cols(dim2) %>% # dim2を結合
rename(Dim2 = `NA`) # カラム名をDim2に変更（デフォルトはNA）
# 結果を出力
head(material2) %>%
kable()
| 石材 | Dim1 | Dim2 |
|---|---|---|
| Obs | 1.8126154 | -0.9823849 |
| Sh | 2.2755113 | -0.3011457 |
| Sh | 2.4030403 | -1.0453552 |
| Sh | 2.0050043 | -1.0033981 |
| Sh | 1.6181339 | -1.3398001 |
| Sh | -0.8991812 | -0.6126366 |
OK.目的達成です。
このデータで再度主成分得点の散布図を出力します。
library(ggthemes)
material2 %>%
ggplot(aes(x = Dim1 ,
y = Dim2 ,
colour = 石材)) +
geom_point( ) +
scale_colour_ptol() +
theme_minimal()
むー、石材ごとに違いが見えるかと思いましたが、そううまくはいきませんでしたね。
頁岩（Sh）は、幅広タイプとノーマルタイプの2種類があるってことですかね？
以上、主成分分析でした。