Cluster_Tutorial

3964 days ago by minami

準備

SageでRを使用するために必要な前準備を以下で行います。
# Rのユーティリティ関数を読み込む attach(DATA+'RUtil.py') 
       

クラスタリング

データがどのようなグループに分かれているか調べる方法をクラスタリングと呼びます。

以下の手順でクラスタリングの手法について整理します。

  1. 準備として、Excelシートからのデータ読み込み
  2. ユークリッド距離を使った階層型クラスタリング
  3. k-means法によるクラスタリング
  4. クラスリング結果を2次元で表示
  5. 最適なクラスターの数の求め方
  6. クラスタリング結果を元データの分布で見る

Excelシートからのデータ読み込み

RでExcelシートからデータを読み込むには、gdataパッケージのread.xls関数を使用します。
# r('install.packages("gdata")') r('library(gdata)') 
        
[1] "gdata"     "stats"     "graphics"  "grDevices" "utils"    
"datasets"  "methods"   "base"     
[1] "gdata"     "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"

ここでは、ノートブックにある成績のデータ（seiseki.xls）を読み込みます。

row.names=1の引数で、行の名称を1列目の値とすることを指定します。
filename = DATA + 'seiseki.xls' r('seiseki <- read.xls("%s", sheet=1, row.names=1)'%filename) r('seiseki') 
        
     算数 理科 国語 英語 社会
田中   89   90   67   46   50
佐藤   57   70   80   85   90
鈴木   80   90   35   40   50
本田   40   60   50   45   55
川端   78   85   45   55   60
吉野   55   65   80   75   85
斎藤   90   85   88   92   95
     算数 理科 国語 英語 社会
田中   89   90   67   46   50
佐藤   57   70   80   85   90
鈴木   80   90   35   40   50
本田   40   60   50   45   55
川端   78   85   45   55   60
吉野   55   65   80   75   85
斎藤   90   85   88   92   95

距離を計算する

クラスタリングするために、教科毎の点数の差の自乗和から生徒の距離関係を求めます。

距離が小さいほど類似しており、距離が大きいほど異なっていると判断します。
# ユークリッド距離を計算 r('seiseki.d <- dist(seiseki)') # 結果を整数に丸めて表示 r('round(seiseki.d)') 
        
     田中 佐藤 鈴木 本田 川端 吉野
佐藤   69                         
鈴木   34   81                    
本田   60   64   53               
川端   28   61   21   47          
吉野   63   12   76   54   56     
斎藤   68   38   88   92   68   46
     田中 佐藤 鈴木 本田 川端 吉野
佐藤   69                         
鈴木   34   81                    
本田   60   64   53               
川端   28   61   21   47          
吉野   63   12   76   54   56     
斎藤   68   38   88   92   68   46

階層クラスタリング

階層クラスタリングは、hclust関数を使って計算します。
r('sei.hc <- hclust(seiseki.d)') 
        
Call:
hclust(d = seiseki.d)

Cluster method   : complete 
Distance         : euclidean 
Number of objects: 7 
Call:
hclust(d = seiseki.d)

Cluster method   : complete 
Distance         : euclidean 
Number of objects: 7

計算結果をデンドログラムにプロットします。

類似しているものが近くに集まってプロットされます。

残念ながら、Mac OS Xでは日本語が文字化けします。
graph = preGraph("fig5_4.pdf") r('plot(sei.hc)') postGraph(graph, fac=0.7) 
       

氏名を数字の1から7に置き換えてプロットし直しました。

大きく分けて2つのグループに分かれているのが、見て取れます。
# 人の名前を数字に置換して再度グラフを表示 print r('rownames(seiseki)') r('rownames(seiseki) <- c(1:7)') print r('rownames(seiseki)') r('seiseki.d <- dist(seiseki)') r('sei.hc <- hclust(seiseki.d)') graph = preGraph("fig5_4a.pdf") r('plot(sei.hc)') postGraph(graph, fac=0.7) 
        
[1] "田中" "佐藤" "鈴木" "本田" "川端" "吉野" "斎藤"
[1] "1" "2" "3" "4" "5" "6" "7"

[1] "田中" "佐藤" "鈴木" "本田" "川端" "吉野" "斎藤"
[1] "1" "2" "3" "4" "5" "6" "7"

k-means法

階層クラスタリングは、直感的とよくマッチするので分かりやすいのですが、要素が多くなると処理に時間が掛かります。

次にK個のグループ分けを逐次行うk-means法を使って、クラスタリングをします。

k-means法のクラスタリングには、kmeans関数を使用します。

seiseki.km$clusterにどのクラスに割り振られたかを示す番号がセットされます。

佐藤（2）、吉野（6）、斎藤（7）のグループ１と

田中（1）、鈴木（3）、本田（4）、川端（5）のグループ２に分けられ、階層クラスタリングの結果と一致します。

 
# k-meansでクラスタリング r('seiseki.km <- kmeans(seiseki,2)') r('seiseki.km$cluster') 
        
1 2 3 4 5 6 7 
2 1 2 2 2 1 1 
1 2 3 4 5 6 7 
2 1 2 2 2 1 1

クラスタリング結果の可視化

多次元の属性から求められたグループ分けがどのような関係になっているかを2次元にプロットするのが、clusplot関数です。

これによってグループ内のメンバーの相対的な位置づけやグループ間の距離が視覚的に理解できます。
# 分布のかたまりを2次元のグラフで表示する r('library(cluster)') graph = preGraph("seiseki_cluster.pdf") r('clusplot(seiseki, seiseki.km$cluster)') postGraph(graph, fac=0.7) 
       

クラスタリングの結果と属性内での分布をPairsグラフで表示するのが、clPairs関数です。

成績の結果をみると、2つのグループは理系（赤）と文系（青）であることがMathematicsとScienceの分布で分かります。
# 元のデータでのクラスタリングが行われているか表示する # r('install.packages("mclust")') r('library(mclust)') # 科目名を以下の様に変更 r('colnames(seiseki) <- c("Mathematics", "Science", "Japanese", "English", "Sociery")') graph = preGraph("seiseki_clPairs.pdf") r('clPairs(seiseki, seiseki.km$cluster)') postGraph(graph, fac=0.7) 
       

クラスターの数を求める

k-means法は、とても便利で高速な手法ですが、肝心のクラスターの数をいくつにするのが良いのか、その指標が必要になります。

そこで、混合分布によるクラスター分析手法を使って最適なクラスター数を求めます。

mclustパッケージのmclustBIC関数を使って最適なクラスター数を計算したのが、以下の図です。

BICの値が高いほど良いモデルの指標です。以下の例ではVEVモデル、VVVモデルが最もよくクラスターの数も２、３が良いことが求まります。
# k-meansのクラスター数をどのように求めるのか # irisのデータを元に計算 graph = preGraph("5_8.pdf") r('irisBIC <- mclustBIC(iris[,1:4])') r('plot(irisBIC)') postGraph(graph, fac=0.7) 
       

irisのデータは品種の異なるユリ（setosa, versicolor, virginica）の花の萼(がく)の長さと幅、花弁の長さと幅を測ったもので、それぞれ50個のサンプルが入っています。

データの最後に品種名がセットされています。

hc関数で混合分布モデルのクラスタリングを行い、hclass関数でクラス分けをします。

結果として、setosa, versicolorは、100%分類できていますが、virginicaは、50個中14が誤識別されています。
# クラスター数2 or 3で一番BICが大きくなっているので、k=3でクラス分けしてみる # VEVは実装されていないので、VVVを使用する # 混合分布VVVモデルを使って階層クラスタリングを行う r('iris.hc <- hc(modelName = "VVV", data = iris[,1:4])') # クラス数を3でクラス分けする r('iris.hcl <- hclass(iris.hc, 3)') # クラス分けされた結果と正解のクロス集計を取り、virginicaが14個、誤識別されている r('table(iris[,5], iris.hcl)') 
        
            iris.hcl
              1  2  3
  setosa     50  0  0
  versicolor  0 50  0
  virginica   0 14 36
            iris.hcl
              1  2  3
  setosa     50  0  0
  versicolor  0 50  0
  virginica   0 14 36

<h3>irisクラスタリング結果の可視化</h3>

irisデータの可視化の結果を以下に示します。

Sepal.Length, Sepal.Width属性においてVirginica, versicolorの識別が混ざっていることが見て取れます。
# graph = preGraph("5_9.pdf") r('irisBIC <- mclustBIC(iris[,1:4])') r('clPairs(iris[,1:4], cl=iris.hcl)') postGraph(graph, fac=0.7) 
       
graph = preGraph("5_10.pdf") r('clusplot(iris[,1:4], iris.hcl)') postGraph(graph, fac=0.7)