Rのhclust()によるデンドログラムの作成（ウォード法）
最終更新日 ： 2016.9.17
クラスター分析との応用
クラスタリング（clustering）とは「分類する」という意味であるが、少し表現を変えると「ある個体が特定のクラスターに所属するかどうかを判定する」ということであり、そのための統計解析手法がクラスター分析（cluster analysis）である。これまで説明してきた分散分析モデル、回帰モデル、主成分分析モデル、因子分析モデルなどはいずれにせよ Y = X1 + X2 + ... + Xp といった線形モデルとして表現することができた。これらは「ある説明変数によって目的変数を説明する（または予測する）」ということを表現していた。しかしクラスター分析は個体を分類するという観点での分析手法であり、いままでの手法とは目的が異なるものといえる。
なお冒頭の説明で「分類する」という表現と「判定する」という表現を用いたが、判定は判別と言い換えてもよく、すでに説明した判別分析とクラスター分析は手法上まったく異なるが、実際の運用に際しては非常に密接な関係がある。そのため、本ページでは最初にクラスター分析（個体を分類する）を紹介し、すでに紹介した因子分析との組み合わせによる応用的な解析法、そして「未知の個体が」どのクラスターに所属するかを予測するための方法を説明していく。
■クラスター分析 〜個体を分類する〜
ここでは再び、因子分析の解説で用いたデータセット（豊田秀樹「共分散構造分析 Amos編」東京図書（2007）-p64：東京図書のウェブサイトの書籍紹介のページよりダウンロードできるc03test.xls／.datファイルはこちら）を用いることにする。このデータセットは100人分の6教科のテスト得点（国語、地理、英語、数学、化学、物理）が観測されたものである。
この6教科のテスト得点のデータに基づき、学生100人をいくつかのクラスターに分類することを考えてみよう。少し想像力を働かせてみると以下のような仮説が思いつくのではないだろうか。
クラスター分析において、最初に考えるべきこと（思いつく疑問点）は「いったいいくつのグループにわけられるのだろうか」ということである。適切なクラスター数を決定する方法として「x-means」や「モデルに基づいたクラスター分析」などと呼ばれる手法が存在するが、必ずしも実質的な観点からみて適切なクラスター数を教えてくれるものではない。そのため、適切なクラスター数を決めるには、クラスター分析の結果を解釈することで「実質的に妥当だと考えられるクラスター数」を解析者自身が決めていくしかないというのが実情であろう。
とはいえ、まったく何の情報もないところからクラスター数を決めるというのも難しい。そこで方法の1つとして、階層的クラスター分析によりデンドログラムを作成して、それを参考に大まかなクラスター数を決めるというものを紹介する。Rで階層的クラスター分析を実行するためにhclust()という関数が用意されている。このhclust()には距離行列とクラスタリングの方法を指定することができる。
距離の定義には様々なものがあるが、一般的にユークリッド距離を用いることが多い。これはdist()によって得られる。クラスタリングの方法として、これもまた様々なものがあるが、経験的にウォード法と呼ばれるものがよい分類をしてくれることが知られている（ただしウォード法を用いる際には、ユークリッド距離の自乗を指定する必要がある）。
## データセットの読み込み dat <- read.table("factanal.dat", header=TRUE, row.names=1, sep="\t") d <- dist(dat) # 距離行列の作成 result <- hclust(d^2, method="ward") # ウォード法によるクラスター分析 plot(result) # デンドログラムを描画する
Rのhclust()によるデンドログラムの作成（ウォード法）
さて、このトーナメント表のような図のことをデンドログラムというが、一見すると枝分かれしたその先に葉っぱが生茂っているようにみえる。この生茂った葉っぱのかたまりが「クラスター」と呼ばれるものである。階層的クラスター分析では、この枝を「ある高さで切る」ことにより、その切られた枝に含まれる葉っぱが1つのクラスターになる（下図を参照）。
ある高さで「枝」を切ったときのクラスター（上図ではクラスター数が4である）
どの高さで枝を切るかは解析者が決めることになるが、どの高さで切るのが適切かはやはり難しい問題であるが、少なくともどのくらいのクラスター数になりそうかを視覚的に（てっとり早く）把握するためには役立つであろう。例えば、これよりも低い位置で枝を切ろうとするとクラスター数は多すぎるかもしれない。逆にこれよりも高い位置で枝をきってしまうと、クラスター数が少なすぎるかもしれない。とりあえず今回はクラスター数が4であることを仮定してみることにする。
ところで、この階層的クラスター分析はデンドログラムによって視覚的に把握できる利点があるものの、変数の数（次元数）とサンプルサイズ（個体の数）が多くなると計算負荷が高くなり膨大な時間がかかるという問題がある。そこで高速でかつ経験的によい分類をしてくれるとされている方法にk-means法と呼ばれる手法がある。このk-means法はクラスター分析の中でも特に非階層的クラスター分析と呼ばれている。
また（階層的クラスター分析を行うhclust()でもできなくはないが）Rに実装されているkmeans()という関数を用いてk-means法によるクラスター分析を行うことで、クラスターの解釈を行うために有用なアウトプットを簡単に得ることができる。以下にkmeans()の実行例を示す。
## kmeans法によるクラスター分析。クラスター数を「4」として実行 result2 <- kmeans(dat, centers=4, iter.max=200) ## kmeans()が返した結果オブジェクトの中身を確認 names(result2) [1] "cluster" "centers" "totss" "withinss" "tot.withinss" [6] "betweenss" "size" ## centersの中に「クラスター平均」が格納されている result2$centers 国語 地理 英語 数学 化学 物理 1 44.33333 56.00000 52.33333 51.33333 51.33333 43.00000 2 59.72973 64.45946 65.27027 79.86486 73.91892 66.62162 3 63.75000 67.70833 62.70833 66.66667 41.66667 60.41667 4 66.66667 64.79167 64.79167 55.41667 67.91667 55.20833 ## clusterの中に各個体の所属クラスター番号が格納されている result2$cluster 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 4 2 2 3 1 4 2 3 1 3 3 2 1 2 4 3 2 3 2 2 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 1 3 3 1 1 2 2 2 2 4 2 4 2 1 4 2 2 3 2 4 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 2 2 3 1 3 3 4 3 3 4 4 3 2 3 1 4 2 4 1 4 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 3 2 1 3 4 2 2 2 2 4 1 4 3 3 1 2 2 3 2 2 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 2 4 3 4 4 2 2 2 2 4 2 4 4 1 1 4 4 3 2 3
各クラスターがどのような意味を持っているか（どのような特徴を持ったグループであるあるのか）はクラスター平均を読み解くことで把握することができる。上記の result2$centers に格納されているクラスター平均を分かりやすく表にまとめて色塗りをしてみよう（下表を参照）。クラスター平均表とは、各クラスターにおける科目ごとの平均得点をまとめたものである。例えばクラスター「1」に所属する人数は15人であり、その15人の国語の得点の平均値が44.3である。また、科目内平均とは各科目のクラスター平均の平均であり、国語の科目内平均は (44.3 + 59.7 + 63.8 + 66.7) / 4 = 58.6 である。一方でクラスター内平均とは各クラスターの科目の平均の平均値であり、クラスター1のクラスター内平均は (44.3 + 56.0 + 52.3 + 51.3 + 51.3 + 43.0) / 6 = 49.7 である。
| クラスター | 所属人数 | 国語 | 地理 | 英語 | 数学 | 化学 | 物理 | クラスター内平均 |
| 1 | 15 | 44.3 | 56.0 | 52.3 | 51.3 | 51.3 | 43.0 | 49.7 |
| 2 | 37 | 59.7 | 64.5 | 65.3 | 79.9 | 73.9 | 66.6 | 68.3 |
| 3 | 24 | 63.8 | 67.7 | 62.7 | 66.7 | 41.7 | 60.4 | 60.5 |
| 4 | 24 | 66.7 | 64.8 | 64.8 | 55.4 | 67.9 | 55.2 | 62.5 |
| 科目内平均 | 58.6 | 63.2 | 61.3 | 63.3 | 58.7 | 56.3 | - |
クラスター数が4であるときのクラスター平均表を見てみると、クラスター1の特徴はどの科目もイマイチ点が低い「落第候補生クラスター」などとラベル付けすることができるかもしれない。クラスター2については、どの科目も平均得点が高いが、特に数学、化学、物理といった科目の平均値が高いことから「理系得意クラスター」と名付けられるかもしれない。ところがクラスター3とクラスター4はどちらも似たような傾向をもっており、いまひとつ解釈がしにくいといえる。このような場合、不要にクラスター数を多く設定していることが考えられる。ためしにクラスター数を「3」にして、再度クラスター分析を行ってみよう。
| クラスター | 所属人数 | 国語 | 地理 | 英語 | 数学 | 化学 | 物理 | クラスター内平均 |
| 1 | 34 | 55.3 | 59.7 | 59.1 | 51.2 | 59.1 | 50.9 | 55.9 |
| 2 | 21 | 65.5 | 68.3 | 64.0 | 68.1 | 40.5 | 57.1 | 60.6 |
| 3 | 45 | 61.1 | 65.3 | 64.6 | 77.4 | 72.8 | 65.7 | 67.8 |
| 科目内平均 | 60.6 | 64.5 | 62.6 | 65.6 | 57.5 | 57.9 | - |
クラスター1の特徴は先ほどの結果と同じく「落第候補生クラスター」のようである。クラスター2は国語、地理、英語の文系科目の平均値が高いため「文系得意クラスター」といえるかもしれない。こちらは先ほどのクラスター3・4の特徴を表現しているといえそうである。そしてクラスター3は「理系得意クラスター」といえそうである。
実際には、より注意深く見てみると、文系得意クラスター（クラスター2）であっても数学の得点は科目内平均の65.6よりも高くなっていたり、理系得意クラスター（クラスター3）における英語の平均が科目内平均よりも高くなっていたりする。クラスターの解釈で難しいのは、単に「科目内平均よりも高い」ということが重要であるとは限らず、クラスター内平均、科目内平均、クラスター×科目の個別の平均値との比較など総合的に解釈することにある。そういう意味では「理系の得意な学生は一般的に英語も得意である」ということを表しているかもしれないのである。その背景には「科学論文は英語で読み書きしなければならない」という事情があるのかもしれない。こうした実質的な観点も含めて解釈をすることにクラスター分析の難しさが存在するといえる（わるい言い方かもしれないが、クラスター分析は「結果オーライ」な手法であるともいえる）。
■因子得点を利用したクラスター分析
クラスター分析を行う際に難しいのは解析すること自体というよりは、むしろ解析結果から各クラスターへの意味づけを行うことにあるということは分かってもらえただろうか。上述したような例は非常に単純な問題であるが、より多くの変数に基づいてクラスタリングを行う場合にはその解釈がより困難になる。解析上の観点でも「次元の呪い」とよばれる問題が生じるが、そもそもクラスター解析で重要なことは意味のあるクラスターを作成することにあるので、解釈可能な変数を用いる（解釈可能な変数の数を考える）べきである。一般に機械学習の分野では、単に次元を圧縮するという目的（非常に解析的な観点に偏った考え方である）で主成分分析を行い、その主成分得点を用いてクラスタリングを行う方法などが提案されている。しかしながら、前述したように主成分分析のような次元圧縮はあくまでも「解釈を容易にする」という目的で実行されるべきである。
ここに主成分分析の例をあげたが、今回は因子分析との組み合わせによるクラスター分析の応用事例を示すことにする。まずは因子分析のページで説明したように、先ほど用いた6科目（100人分）のデータに対して2つの因子得点を算出する。その因子得点（第1因子：文系科目、第2因子：理系科目の2変数（2次元）のデータ）に対してk-means法によるクラスター分析を行ってみる。
## 因子分析を実行する fa.result <- factanal(dat, factors=2, rotation="promax", scores="regression") ## 因子負荷量行列を表示 fa.result Loadings: Factor1 Factor2 国語 0.470 地理 -0.116 0.614 英語 0.555 数学 0.523 0.104 化学 0.652 -0.209 物理 0.536 0.128 ## factanal()の実行結果オブジェクトの中身 names(fa.result) [1] "converged" "loadings" "uniquenesses" "correlation" "criteria" [6] "factors" "dof" "method" "rotmat" "scores" [11] "STATISTIC" "PVAL" "n.obs" "call" ## 因子得点はscoresに格納されている fs <- fa.result$scores ## 因子得点の最初の10行を確認 fs[1:10, ] Factor1 Factor2 1 -0.2574858 0.3233926 2 0.4760993 -0.4694338 3 0.3334185 0.7696370 4 -1.1862578 1.5904724 5 -1.0501571 -0.7166635 6 -0.2216984 0.3289725 7 0.8613351 0.3009074 8 -1.7100095 1.4013709 9 -1.2888163 -0.9058265 10 -1.4974776 0.4650623 ## 因子得点を用いてクラスタリングを行う km.result <- kmeans(fs, 3, 200) ## クラスター平均 km.result$centers Factor1 Factor2 1 -0.4438874 0.8546890 2 -0.4793381 -0.6485683 3 0.8176907 -0.3084898 ## 各クラスターの所属人数 table(km.result$cluster) 1 2 3 35 29 36
この結果を、今までと同じようにクラスター平均表にまとめてみよう（下表）。これを看ると6科目の時よりも解釈がしやすくなっており、クラスター1が文系得意クラスター、クラスター2が落第候補クラスター、クラスター3が理系得意クラスターであることがわかる。なお因子得点は平均が0.00になるように標準化されているので、この場合は単純に正負の符号だけで判断できてしまう（因子内平均がほぼ0.00なので、プラスならば高得点、マイナスならば低得点ということがわかる）。ただし、解釈がしやすいという利点がある反面、6科目のクラスター分析の結果に比べ情報量が少なくなっている点には留意しなければならない。3つのクラスターの特徴は把握しやすくなったが、例えば「理系の得意な学生は一般的に英語も得意である」といった知見は、この結果からだけでは得ることができない。
| クラスター | 所属人数 | 第1因子：理系科目 | 第2因子：文系科目 | クラスター内平均 |
| 1 | 35 | -0.44 | 0.85 | 0.21 |
| 2 | 29 | -0.48 | -0.65 | -0.56 |
| 3 | 36 | 0.82 | -0.31 | 0.25 |
| 因子内平均 | -0.04 | -0.03 | - |
■未知の個体の所属クラスターを予測する
6科目のデータに基づいて（あるいは2つの因子得点に基づいて）クラスタリングを行った結果、3つのクラスターに分類することができた。ここで例えばNo.101番目の個体（学生）が現れたとき、その個体がどのクラスターに所属するかはどのようにして知ることができるだろうか。このような判別予測をしたい場合に役立つのが既に説明した判別分析である。しかし判別分析のページで説明した問題は、判別するものが男性（male）か女性（female）かという2水準のデータであった。2値の判別の場合にはフィッシャーの線形判別関数を用いた方法や2値ロジスティック回帰モデルを用いることができた。ところが今回は3つのクラスター、つまり3水準のデータを判別予測するという問題である。
3水準以上の判別予測を行う場合にふさわしい統計モデルとして、古典的には多項ロジスティック回帰分析をあげることができるであろう。最近では（といっても、それなりに古くから知られている手法であるが）サポートベクターマシン（support vector machine, SVM）という手法がある。この手法についての詳細は省くが、Rのe1071パッケージにはそれを行うための関数svm()が実装されている。
## 6科目のデータに、kmeansによって作られたクラスター番号を列結合する dat2 <- cbind(dat, km.result$cluster) colnames(dat2)[7] <- "CLS" dat2[1:10, ] 国語 地理 英語 数学 化学 物理 CLS 1 65 65 70 50 65 60 1 2 50 45 80 80 65 55 3 3 60 75 80 85 75 50 3 4 65 85 75 60 40 55 2 5 50 55 55 45 50 40 1 6 70 60 70 50 60 65 1 7 55 75 65 85 70 75 3 8 80 60 80 60 20 45 2 9 30 65 55 55 55 20 1 10 65 65 60 60 30 40 2 ## パッケージe1071にあるsvm()を実行する install.packages("e1071") library(e1071) svm.result <- svm(as.factor(dat2$CLS) ~ ., data=dat2, probability=TRUE, kernel="radial") predict(svm.result, newdata=dat2, probability=TRUE) # 以下は出力結果：svmによる判別予測結果 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 3 3 2 1 1 3 2 1 2 2 3 1 3 1 3 3 2 3 3 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 1 2 2 1 1 3 3 3 3 1 3 3 3 1 1 3 3 2 3 3 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 3 3 2 1 1 2 1 2 2 3 1 2 3 2 1 1 3 2 1 1 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 2 3 1 3 1 3 3 3 3 1 1 1 2 1 1 3 3 2 3 3 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 3 3 2 3 1 3 3 3 3 3 3 1 1 1 1 1 1 2 3 2 # 各クラスターの所属確率が出力される attr(,"probabilities") 1 3 2 1 9.891200e-01 0.0078417869 0.0030382421 2 2.802544e-02 0.9514702523 0.0205043037 3 1.570185e-02 0.9771649530 0.0071331948 4 1.325426e-02 0.0233780187 0.9633677177 5 9.984166e-01 0.0004408377 0.0011425597 6 9.607000e-01 0.0221160237 0.0171840103 (以下省略) ## どのくらい正確に判別予測されているかどうかを確認してみる prd <- predict(svm.result, newdata=dat2) # prdに判別予測結果を格納 table(prd, dat2$CLS) # 実測値と予測値をクロス集計する prd 1 2 3 # 完璧に判別予測できている 1 34 0 0 2 0 21 0 3 0 0 45 ## 未知の個体No.101のデータフレームを作成 no101 <- data.frame(国語=23, 地理=40, 英語=45, 数学=69, 化学=90, 物理=78) no101 国語 地理 英語 数学 化学 物理 1 23 40 45 69 90 78 ## No.101がどのクラスターに所属するかを予測 no101.prd <- predict(svm.result, no101) no101.prd 1 3 Levels: 1 2 3
参考文献
[1] 水野欽司「多変量データ解析広義」朝倉書店（2002）
[2] B. エヴェリット著・石田基広ほか訳「RとS-PLUSによる多変量解析」シュプリンガー・ジャパン（2007）
[3] 豊田秀樹「データマイニング入門」東京図書（2009）
[4] 豊田秀樹「共分散構造分析 Amos編」東京図書（2007）
[5] 阿部重夫「パターン認識のためのサポートベクトルマシン入門」森北出版（2011）
