多変量データの扱い
高校生のための統計学習教材（総務省統計局）
の解説です。
上記サイトの第２講「多変量データの扱い」をダウンロードし、テキストとしてください。
テキスト（p9〜,p19〜）
にそって統計ソフト「RStudio」の適用例を提示します。必要に応じて「なるほど統計学園」の参考ページを参照してください。
1 主成分分析
物事にはたくさんの要因（変数）があります。しかし、変数が多いと理解が難しくなります。
そこで、理解しやすいように、多くの変数をまとめることが必要になります。
対象の特徴をデータで表す際、単一の変数だけではなく多くの変数を使うことで対象の特徴を詳しく分析することができます。
単回帰分析に続いて重回帰分析が考えられたのはこのためです。
実際、変数の数を増やせば決定係数も大きくなり、当てはまりがよくなるように見えます。
しかし，変数の数が増えすぎると、かえって対象の特徴を総合的に捉えることが難しくなります。
主成分分析は多くの変数から、対象を最もよく表現する新しい変数を作成し、対象の特徴を総合的にとらえる分析方法です。
例えば、５教科のテスト成績は５次元のデータですが、５教科の点数をもとに「理系学力」「文系学力」という変数を作ればデータは２次元になり、学力の特徴を把握しやすくなります。
特徴は分散が大きいほど識別しやすくなります。対象が離れて見える方向から観察することが重要です。
このため主成分となる変数は、
「新変数の分散が大きくなる」「新変数間の相関がなく情報が重複しない」
ようにつくられます。
次に示す例は、47都道府県
31品目の食料品消費傾向の分析です。（総務省統計局）
各都道府県の特徴はどうなっているでしょうか？
31次元のデータの全体像を把握することは難しいですが、データを２次元に圧縮し図示できれば、
各都道府県の特徴を視覚化できます。
まずデータファイルを読み込み、データフレーム化します。
読み込んだ項目は次の通りです。
## [1] "都道府県" "米【円】"
## [3] "パン【円】" "麺類【円】"
## [5] "他の穀類【円】" "生鮮魚介【円】"
## [7] "塩干魚介【円】" "魚肉練製品【円】"
## [9] "他の魚介加工品【円】" "生鮮肉【円】"
## [11] "加工肉【円】" "牛乳【円】"
## [13] "乳製品【円】" "卵【円】"
## [15] "生鮮野菜【円】" "乾物・海藻【円】"
## [17] "大豆加工品【円】" "他の野菜・海藻加工品【円】"
## [19] "生鮮果物【円】" "果物加工品【円】"
## [21] "油脂【円】" "調味料【円】"
## [23] "菓子類【円】" "主食的調理食品【円】"
## [25] "他の調理食品【円】" "茶類【円】"
## [27] "コーヒー・ココア【円】" "他の飲料【円】"
## [29] "酒類【円】" "一般外食【円】"
## [31] "学校給食【円】"都道府県名を行タイトルにします。
30品目について、各都道府県の支出額は次の通りです。
全体像をつかむため、品目ごとの箱ヒゲ図を作成します。
## 'data.frame': 47 obs. of 31 variables:
## $ 都道府県 : chr "北海道" "青森県" "岩手県" "宮城県" ...
## $ 米【円】 : num 2113 1742 2004 2090 2550 ...
## $ パン【円】 : num 1760 1617 1633 1632 1296 ...
## $ 麺類【円】 : num 1017 1163 1111 1114 1173 ...
## $ 他の穀類【円】 : num 258 204 264 213 152 304 235 280 315 272 ...
## $ 生鮮魚介【円】 : num 3073 3379 3447 3124 3423 ...
## $ 塩干魚介【円】 : num 923 1056 1022 902 1189 ...
## $ 魚肉練製品【円】 : num 434 371 474 563 385 468 465 415 470 470 ...
## $ 他の魚介加工品【円】 : num 707 737 709 634 831 704 692 676 759 806 ...
## $ 生鮮肉【円】 : num 3468 3513 3358 3454 3721 ...
## $ 加工肉【円】 : num 1075 1051 929 999 1018 ...
## $ 牛乳【円】 : num 915 922 1112 1112 998 ...
## $ 乳製品【円】 : num 1087 990 1267 1213 976 ...
## $ 卵【円】 : num 477 461 581 564 541 610 595 562 565 567 ...
## $ 生鮮野菜【円】 : num 4441 4465 4795 4959 4786 ...
## $ 乾物・海藻【円】 : num 426 496 533 520 530 551 454 486 509 505 ...
## $ 大豆加工品【円】 : num 734 828 1103 1011 890 ...
## $ 他の野菜・海藻加工品【円】: num 800 783 811 889 840 933 889 871 907 873 ...
## $ 生鮮果物【円】 : num 2157 2054 2514 2575 2362 ...
## $ 果物加工品【円】 : num 120 119 91 134 90 122 115 125 150 140 ...
## $ 油脂【円】 : num 258 238 216 211 219 236 227 249 277 242 ...
## $ 調味料【円】 : num 2402 2320 2476 2467 2477 ...
## $ 菓子類【円】 : num 3980 3948 4397 4490 4155 ...
## $ 主食的調理食品【円】 : num 2301 2783 2635 3265 2708 ...
## $ 他の調理食品【円】 : num 3331 4165 4156 4264 4095 ...
## $ 茶類【円】 : num 704 626 730 798 758 ...
## $ コーヒー・ココア【円】 : num 905 1166 891 817 826 ...
## $ 他の飲料【円】 : num 1744 2172 1983 1803 1881 ...
## $ 酒類【円】 : num 3148 4126 3504 2800 3830 ...
## $ 一般外食【円】 : num 8834 6881 8104 9724 9634 ...
## $ 学校給食【円】 : num 688 804 665 710 814 878 969 954 916 891 ...
最後から２番目の「一般外食」は金額が桁違いに大きいいことが分かります。
主成分分析を行う関数は “prcomp”です。
３０次元の支出を２次元に圧縮し、biplot関数で図示することが可能になります、
# 主成分分析
par(family = "HiraKakuProN-W3") #日本語フォントの指定
par(ps = 6)
result<-prcomp(x1,scale=F)
biplot(result) 
外食に比べると、単食料品は金額が小さいため、(0,0)付近に密集して見えます。
赤い文字は、その品目についての第１主成分PC1と第２主成分PC２の構成比を、点としてプロットしたものです。
東京、神奈川は、主成分の比率が似通っており、外食影響が大きいことが分かります。
福井、富山は支出傾向が似通っており、外食も少ないことが読み取れます。
par(family = "HiraKakuProN-W3")
plot(result$x[,1], result$x[,2], pch=19, cex=0.6, xlab="PC1", ylab="PC2")
text(result$x[,1], result$x[,2], rownames(x1), cex=0.6, pos=4)
主成分分析では、多次元データを２次元に圧縮して図示することにより、対象の類似性などを把握することができます。
2 クラスター分析
コロナ対策でクラスターという言葉が有名になりました。もとはブドウの房のようにサンプルが集まったものを言います。パンデミックの初期においては、感染者がブドウの房のように固まって現れ、その房を摘み取ることが効果的な感染対策になった訳です。
データサイエンスでいうクラスターは類似したサンプルをまとめた塊を言います。サンプルをクラスターに分類することは、全体の傾向の把握に役立ちます。
類似度として最も基本的なものはユークリッド距離です。つまり、距離の近いものをまとめていくことになります。
2.1 階層的クラスタリング
全国の都道府県を消費傾向によってクラスターに分けてみましょう。 都道府県の数を8つ、食品項目4つにして類似した都道府県をまとめます。
#クラスター分析
(x2<- dt1[c(1:9),c(2:5)])#都道府県と食品項目を絞る| 米【円】 | パン【円】 | 麺類【円】 | 他の穀類【円】 | |
|---|---|---|---|---|
| 北海道 | 2113 | 1760 | 1017 | 258 |
| 青森県 | 1742 | 1617 | 1163 | 204 |
| 岩手県 | 2004 | 1633 | 1111 | 264 |
| 宮城県 | 2090 | 1632 | 1114 | 213 |
| 秋田県 | 2550 | 1296 | 1173 | 152 |
| 山形県 | 2443 | 1627 | 1352 | 304 |
| 福島県 | 2528 | 1619 | 1167 | 235 |
| 茨城県 | 2237 | 1781 | 1168 | 280 |
| 栃木県 | 2146 | 1779 | 1277 | 315 |
(dat.d <- dist(x2))#都道府県同士のユークリッド距離を求め行列で表す## 北海道 青森県 岩手県 宮城県 秋田県 山形県 福島県 茨城県
## 青森県 426.9918
## 岩手県 192.0469 274.2335
## 宮城県 168.3657 351.8679 100.0350
## 秋田県 664.7082 871.0390 654.2729 575.9323
## 山形県 490.8462 732.9543 502.4321 435.3837 419.7082
## 福島県 463.8265 786.6238 527.9669 441.9344 334.2723 215.1162
## 茨城県 197.7423 526.9934 282.3083 226.3073 591.2724 317.1498 336.0818
## 栃木県 268.8847 463.4404 267.6509 248.4311 658.7033 342.1389 435.9174 146.2566北海道と青森県の距離は
\(\sqrt{(2113-1742)^2+(1760-1617)^2+(1017-1163)^2+(258-204)^2}=426.9918\)
と計算されます。
距離の最も近いものから一つずつまとめてクラスターを作ります。
クラスター間の距離の計算には様々な方法がありますが、省略した場合は”complete
linkage method”が使用されます。ここでは詳細には立ち入りません。
クラスタリングを実行するには、距離行列に関数 hclust を適用します。
res <- hclust(dat.d)
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(res,hang = -1)
縦軸は、クラスター間の距離を表します。
この図はデンドログラムと言われます。デンドログラムをもとに2つのグループに分けると、
par(family = "HiraKakuProN-W3")
plot(res, hang = -1)
rect.hclust(res, k = 2)
デンドログラムを横に描くこともできます。
par(family = "HiraKakuProN-W3")
plot(as.dendrogram(res), horiz = TRUE)#横
47都道府県について、全食品項目について同様のクラスタリングを実施します。
par(family = "HiraKakuProN-W3") #日本語フォントの指定
x1<- dt1[,-1]
dat.d <- dist(x1)
res <- hclust(dat.d)
plot(res,hang = -1,cex=0.5)
2.2 非階層的クラスタリング
k-means法の手順は、
- あらかじめクラスター数kを決める。
- 各サンプルが所属する初期クラスターをランダム決める。
- 初期クラスターの重心(center)の座標を求める。
- 各サンプルを最も近いcenterのクラスターに所属させる。
このように初期クラスターを更新していき、更新しても変わらなくなるまで繰り返します。
クラスター数を２としてクラスターを作ると、
#k-means法（クラスター数は2）
dat.km <- kmeans(x2, center = 2)
#k-means法によるクラスタリング結果の確認
dat.km$cluster## 北海道 青森県 岩手県 宮城県 秋田県 山形県 福島県 茨城県 栃木県
## 2 2 2 2 1 1 1 2 2#パッケージの読み込み
library("cluster")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
#k-means法によるクラスタリング結果の可視化
#clusplot(x2, dat.km$cluster)
#k-means法によるクラスタリング結果の可視化
clusplot(x2, dat.km$cluster, labels = 2)
47都道府県について、6つのクラスターに分けると、
#k-means法（クラスター数は6）
dat.km <- kmeans(x1, center = 6)
#k-means法によるクラスタリング結果の確認
#dat.km$cluster
#パッケージの読み込み
library("cluster")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
#k-means法によるクラスタリング結果の可視化
#clusplot(x1, dat.km$cluster)
#k-means法によるクラスタリング結果の可視化
clusplot(x1, dat.km$cluster, labels = 2,cex=0.6)
初期クラスターをランダムに決めるので、実行するたびに結果が同じとは限りません。
k-means法の初期値依存と言います。
R言語でクラスタリングしてみた を参照させていただきました。