数式なしのクラスター分析－Rでの最適なクラスタ数の予測から分析までの手順

何らかの形でデータを分類しなければならないとしよう。例えばドラクエのようなRPGに登場するモンスターをタイプごとに分類するならば、

• HP、MP
• ゴールド
• 経験値
• つよさ、すばやさ

などなど。どのような項目、特徴で分類されるのか、最終的に何種類に分類されるのか定かではないデータを、何らかの形で分類する方法の一つがクラスター分析だ。

クラスター分析は機械学習で言うところの、教師なし学習に相当する。何らかの特徴を見出し、その特徴に基づいて分類していく。何通りのクラスターに分類されるのかは予想できないが、見当をつけることはできる。

この投稿では最適なクラスター数の特定から、クラスター分析結果の検討まで、Rを用いた一連の手順を紹介する。分析対象データはRに含まれているものを用いるので、RStudioだけで手順を再現することができる。ここでは次の環境を用いている。

• R 3.6.1
• RStudio 1.2.5001

ちなみに現時点でのR、最新版は3.6.2だ。RStudioと共に参照先は、この投稿末尾にまとめている。

豆知識

ｋ	クラスター分割数
クラスター内変動	クラスター内距離、誤差 クラスタ内の全データから、クラスタの重心までの距離。

エルボー法

クラスター内変動（距離、誤差）の平方和を求める。
ｋが増加するほど、クラスター内変動は下がり続ける。その特性から、ｋとその変動をプロットすると変動が急降下した後、変動が安定するポイントがある。この点が肘（エルボー）に相当する。肘に当たるkを、最適なクラスター数とする。
問題は、肘が明確に描画されるプロットは滅多に現れないこと。

シルエット分析法

凝集度	クラスタ内の全データ同士の平均距離
乖離度	隣接するクラスタ同士の全データの平均距離
シルエット・スコア	（乖離度－凝集度）／最大値（乖離度 or 凝集度）

スコア
1 ~ 0	所属するクラスタの中心に近い。 隣接するクラスタから遠い。	分離性能が良い。
= 0	クラスタ間の決定境界付近に存在する。	クラスタリングの分離性能が悪い。
-1 ~ 0		間違ったクラスタに所属している可能性がある。

平均シルエット・スコアを計算し、最もスコアが良好なｋを、最適なクラスター数とする。
問題は、あくまでも平均値であることから、必ずしも分類結果が良好とは限らないこと。

ギャップ統計法

対象データと、その同範囲で一様分布（クラスタ数１）する架空のデータをクラスタリングし、それぞれのクラスタ同士の変動を比較する（これがギャップ）。そのギャップを最大化する最小のｋを、最適なクラスタ数とする。
プロットの見た目に関わらず、機械的に最適値が出力されるため、最もあてにしやすい。

ライブラリとデータの準備

この投稿では、次のライブラリを用いる。

NbClust	クラスタ分析に用いる。
factoextra	クラスタ分析の視覚化に用いる。
Rmisc	プロットの描画支援に用いる。
LifeCycleSavings	今回の分析対象データ。

LifeCycleSavingsは可処分所得に対する貯蓄率のデータだ。国別に集計されている。変数は次のように定義されている。数値は、1960～1970年の平均だ。
いずれについても、この投稿末尾に参照先をまとめている。

データの読み込みを含め、これらを利用するためのコードは次のようになる。データの内dpi（可処分所得に対する金額）だけが百分率ではないため、scaleで標準化している。

R code

library(NbClust)
library(factoextra)
library(Rmisc)

mydata = scale(LifeCycleSavings)

最適なクラスタ数の特定

クラスタ分析を始める前に、最適なクラスタ数を特定する。階層的、非階層的クラスタリングを問わず、多種多様な指標となり得る類似度（距離）、統計量を一括計算し、多数決の仕組みで最適なクラスタ数を提案してくれる関数が、NbClustだ。その指標にはシルエット分析やエルボー法、ギャップ統計量といった、一般的な指標だけでなく、あまり説明されることのない指標も含まれている。
凝集型階層的クラスタリングではウォード法を用いて、非階層的クラスタリングではＫ平均法を用いて計算をする。それぞれから得た結果をfviz_nbclusterで描画する。

R code

myAHCnum = NbClust(mydata, method = "ward.D", index = "all")
myNHCnum = NbClust(mydata, method = "kmeans", index = "alllong")

fig1 = fviz_nbclust(myAHCnum, method = "silhouette")
fig2 = fviz_nbclust(myNHCnum, method = "gap_stat", nboot = 100)

multiplot(fig1, fig2)

重複した出力は省略する。NbClustの提案する最適なクラスタ数は２だ。
multiplotで多数決の結果を表示している。上が階層的クラスタリングの場合、下が非階層的クラスタリングの場合だ。6分割、15分割の分析も比較検討してみる価値がありそうだ。

よくあるエルボー法、シルエット法、ギャップ統計量のプロットも表示してみる。ここでは２～３が最適であると示唆されている。ここでは３は無視する。

R code

fig3 = fviz_nbclust(mydata, FUNcluster = hcut, method = "silhouette", verbose = FALSE)
fig4 = fviz_nbclust(mydata, FUNcluster = hcut, method = "wss", verbose = FALSE)
fig5 = fviz_nbclust(mydata, FUNcluster = kmeans, method = "gap_stat", nboot = 100, verbose = FALSE)

multiplot(fig3, fig4, fig5)

クラスター分析

以上の検討結果をもとに、クラスタ分析を実施する。ここでは各分析結果で以下３つの図を描画している。

1. 階層的クラスタリングの樹形図
2. 階層的クラスタリングのシルエット図
3. 非階層的クラスタリング（Ｋ平均法）のクラスタ図

各図面間でカラー・パレットは共通なのだが、適用順序が異なるため、各図において同色クラスタが同一クラスタを示しているわけではないことに注意してほしい。

R code

mykmeans = function(c, d){
  myAHC = hcut(d, k = c, stand = TRUE, graph = FALSE)
  myNHC = kmeans(scale(d), c, iter.max = 100, nstart = nrow(d))

  fig_a = fviz_dend(myAHC, rect = TRUE)
  fig_b = fviz_silhouette(myAHC, label = TRUE, rotate = TRUE, print.summary = FALSE)
  fig_c = fviz_cluster(myNHC, data = d)
  
  multiplot(fig_a)
  multiplot(fig_b)
  multiplot(fig_c)
}

for (i in c(2, 6, 15)){
  mykmeans(i, mydata)
}

k = 2の場合

k = 6の場合

k = 15の場合

シルエット図で、一部のデータが負の値を記録しているのが分かる。これは分類が誤っている可能性を示唆している。大変見づらいのだが、クラスタ図を見るとクラスタ同士が接しそうな所に位置しているデータを確認できる。
階層的、非階層的を問わず、２つのクラスタ分割結果が良好であることが分かる。

余談

クラスタ分析のプロットにおける水平軸（Dim1）、垂直軸（Dim2）は、主成分分析により導かれる第1、第2成分だ。次の投稿で、主成分分析を取り扱っている。
また、それに続く投稿ではクラスタ分析、主成分分析を伴う話題を提供している。関心があれば参照してほしい。
impsbl.hatenablog.jp
impsbl.hatenablog.jp
impsbl.hatenablog.jp
impsbl.hatenablog.jp

参照

www.r-project.org
rstudio.com
www.rdocumentation.org
www.rdocumentation.org
www.rdocumentation.org
www.rdocumentation.org
www.rdocumentation.org

コード全文

GitHub Gist - clustering

gist.github.com