【R言語】階層的クラスタリング結果の可視化及び客観的なクラスタ数の決定方法について
教師なし学習の代表的な手法の1つであるクラスタ分析について個人的な論点整理を兼ねてZumel氏とMount氏共著の"Practical Data Science with R"の第8章をまとめてみます（端折ったり本文にないコードを追加したりしてます）。
ちなみに、第8章は以下の出版社サイトからFREEで入手可能です。
Manning | Practical Data Science with R
目次
データ
使用するパッケージ
> library(dplyr) > library(ggplot2) > library(fpc) > library(reshape2)
データの読込み
> protein <- read.table("C:\\Users\\admin\\data\\protein.txt", sep="\t", header=TRUE) > glimpse(protein) Observations: 25 Variables: 10 $ Country <fctr> Albania, Austria, Belgium, Bulgaria, Czechoslovakia, Den... $ RedMeat <dbl> 10.1, 8.9, 13.5, 7.8, 9.7, 10.6, 8.4, 9.5, 18.0, 10.2, 5.... $ WhiteMeat <dbl> 1.4, 14.0, 9.3, 6.0, 11.4, 10.8, 11.6, 4.9, 9.9, 3.0, 12.... $ Eggs <dbl> 0.5, 4.3, 4.1, 1.6, 2.8, 3.7, 3.7, 2.7, 3.3, 2.8, 2.9, 4.... $ Milk <dbl> 8.9, 19.9, 17.5, 8.3, 12.5, 25.0, 11.1, 33.7, 19.5, 17.6,... $ Fish <dbl> 0.2, 2.1, 4.5, 1.2, 2.0, 9.9, 5.4, 5.8, 5.7, 5.9, 0.3, 2.... $ Cereals <dbl> 42.3, 28.0, 26.6, 56.7, 34.3, 21.9, 24.6, 26.3, 28.1, 41.... $ Starch <dbl> 0.6, 3.6, 5.7, 1.1, 5.0, 4.8, 6.5, 5.1, 4.8, 2.2, 4.0, 6.... $ Nuts <dbl> 5.5, 1.3, 2.1, 3.7, 1.1, 0.7, 0.8, 1.0, 2.4, 7.8, 5.4, 1.... $ Fr.Veg <dbl> 1.7, 4.3, 4.0, 4.2, 4.0, 2.4, 3.6, 1.4, 6.5, 6.5, 4.2, 2....
階層的クラスタ分析
はじめに、階層的クラスタ分析によりデンドログラムを作成し、25ヵ国がどのようにグループ分けできるかを見たいと思います。
> vars.to.use <- colnames(protein)[-1] > pmatrix <- scale(protein[,vars.to.use]) #Countryカラムを除くすべてのカラムを正規化 > d <- dist(pmatrix, method="euclidean") # ユークリッド距離を求める > pfit <- hclust(d, method="ward.D2") # ウォード法を利用 > plot(pfit, labels=protein$Country, hang=-1) # デンドログラムの作成 > rect.hclust(pfit, k=5, border="red") # 5つのクラスタに分類
上図より、赤い枠で囲まれた5つのグループに分類できそうだと思われます。これらのグループは含まれている国は地理的に近い位置にあるといえます。なので、地理的に近い地域の国でクラスタを形成しているならば、日々の食事習慣みも似通ってくるのではないかと推察できます。なので、次にグループごとに食品（赤身肉、白身肉、魚、果物・野菜）のたんぱく質摂取量がどのような傾向を示しているかを確認してみます。
> groups <- cutree(pfit, k=5) > print_clusters <- function(labels, k) { for(i in 1:k) { print(paste("cluster", i)) print(protein[labels==i,c("Country","RedMeat","Fish","Fr.Veg")]) } }
> print_clusters(groups, 5) [1] "cluster 1" Country RedMeat Fish Fr.Veg 1 Albania 10.1 0.2 1.7 4 Bulgaria 7.8 1.2 4.2 18 Romania 6.2 1.0 2.8 25 Yugoslavia 4.4 0.6 3.2 [1] "cluster 2" Country RedMeat Fish Fr.Veg 2 Austria 8.9 2.1 4.3 3 Belgium 13.5 4.5 4.0 9 France 18.0 5.7 6.5 12 Ireland 13.9 2.2 2.9 14 Netherlands 9.5 2.5 3.7 21 Switzerland 13.1 2.3 4.9 22 UK 17.4 4.3 3.3 24 W Germany 11.4 3.4 3.8 [1] "cluster 3" Country RedMeat Fish Fr.Veg 5 Czechoslovakia 9.7 2.0 4.0 7 E Germany 8.4 5.4 3.6 11 Hungary 5.3 0.3 4.2 16 Poland 6.9 3.0 6.6 23 USSR 9.3 3.0 2.9 [1] "cluster 4" Country RedMeat Fish Fr.Veg 6 Denmark 10.6 9.9 2.4 8 Finland 9.5 5.8 1.4 15 Norway 9.4 9.7 2.7 20 Sweden 9.9 7.5 2.0 [1] "cluster 5" Country RedMeat Fish Fr.Veg 10 Greece 10.2 5.9 6.5 13 Italy 9.0 3.4 6.7 17 Portugal 6.2 14.2 7.9 19 Spain 7.1 7.0 7.2
各クラスタの食品ごとのたんぱく質摂取量の傾向をみると、以下の傾向が見えてきます。
- Cluster 2 is made of countries with higher-than-average red meat consumption.
- Cluster 4 contains countries with higher-than-average fish consumption but low produce consumption.
- Cluster 5 contains countries with high fish and produce consumption.
主成分分析による可視化
クラスタ分析の結果を2次元平面上に示すため、主成分分析による第1主成分（PC1）を横軸に、第2主成分（PC2）を縦軸にしてクラスタごとの分布を確認してみます。
> princ <- prcomp(pmatrix)# プロマックス法を利用 > nComp <- 2 > project <- predict(princ, newdata=pmatrix)[,1:nComp] # 主成分得点を算出（princ$xと同じ） > project.plus <- cbind(as.data.frame(project),cluster=as.factor(groups),country=protein$Country) # データフレームを作成 # ggplotによる図示 > ggplot(project.plus, aes(x=PC1, y=PC2)) + > geom_point(aes(shape=cluster, colour = cluster)) + > theme_dark() + > geom_text(aes(label=country),hjust=0, vjust=1)
クラスタ1（Romania, Yugoslavia, Bulgaria, Albania）とクラスタ5（Greece, Italy, Portugal, Spain）はうまく分離していますが、それ以外のクラスタは1か所に集中して重なりあっています。
ちなみに、第2主成分までの累積寄与率は以下のとおりです。
> summary(princ) Importance of components: PC1 PC2 PC3 PC4 PC5 PC6 PC7 Standard deviation 2.0016 1.2787 1.0620 0.9771 0.68106 0.57020 0.52116 Proportion of Variance 0.4452 0.1817 0.1253 0.1061 0.05154 0.03613 0.03018 Cumulative Proportion 0.4452 0.6268 0.7521 0.8582 0.90976 0.94589 0.97607 PC8 PC9 Standard deviation 0.34102 0.31482 Proportion of Variance 0.01292 0.01101 Cumulative Proportion 0.98899 1.00000
クラスタの安定度の評価
形成されたクラスターが本当に意味のある集まりかどうかを確認するためにfpcパッケージのclusterboot関数を使用します*1。
> kbest.p<-5 # デンドログラムからクラスタ数を5とする > cboot.hclust <- clusterboot(pmatrix,clustermethod=hclustCBI,method="ward.D2", k=kbest.p) > bootMean.data <- data.frame(cluster = 1:kbest.p, bootMeans = cboot_hclust$bootmean) > round(bootMean_data, digit=2) cluster bootMeans 1 0.78 2 0.82 3 0.68 4 0.89 5 0.74
ここでbootMeansの値（Jaccard係数の平均値）が0.6より小さな値を有するクラスタは不安定と判断され、0.6以上0.75以下の場合はある程度はデータのパターンを示していると判断され、0.85より大きな値を有するクラスタはとても安定している（意味のあるクラスタを形成している）と判断されるらしい。上記の結果をみると、
- クラスタ4（北欧クラスタ:Denmark,Finland,Norway,Sweden）は意味のあるクラスタを形成している（安定度が高い）
- クラスタ4以外もある程度妥当なクラスタを形成している（クラスタ3が少し安定度が低い）
ということが分かります。
デンドログラムからクラスタ数を5としたが、bootMeansの値が0.6を下回るような極端に低い値を有するクラスタが存在しなかったため、とりあえず極端に意味のないクラスタは形成されてはいないことを示す1つの証拠が得られたといえそうです。
クラスタ数の決定方法
clusterboot() 関数はクラスター数を所与とする手法でした。ここまではデンドログラムからクラスタ数を判断しましたが、ここではクラスタ数の決定方法について考えてみます。
1つはクラスタ内でのデータの凝集性に着目した手法で、クラスタ内のデータポイントの平均をセントロイドとし、各データポイントとの距離（誤差）の平方和として算出されるクラスタ内距離二乗和という指標です。この指標の値が大きいほどクラスタ内のデータポイントが疎であるといえ、クラスタ数が大きくなるほど密になるため指標の値は小さくなります。
クラスタ内距離二乗和の計算
# 2つのベクトルx,yの距離の2乗和を計算するsqr_edis関数を作成 sqr_edist <- function(x, y) { sum((x-y)^2) } # クラスタ内でのデータポイントとセントロイドの距離の2乗和（クラスタ内距離二乗和）を計算するwss.cluster関数を作成 wss.cluster <- function(clustermat) { c0 <- apply(clustermat, 2, FUN=mean) # クラスタのセントロイドを計算 sum(apply(clustermat, 1, FUN=function(row){sqr_edist(row,c0)})) } # すべてのクラスタ内距離二乗和を計算するwss.total関数を作成 wss.total <- function(dmatrix, labels) { wsstot <- 0 k <- length(unique(labels)) for(i in 1:k){ wsstot <- wsstot + wss.cluster(subset(dmatrix, labels==i)) } wsstot }
エルボー曲線
縦軸にクラスタ内距離二乗和の値、横軸にクラスタ数をとりプロットしたエルボー曲線を見てみる。クラスタ内距離二乗和はクラスタ数が多くなるほど減少するため、減少幅が一定の幅に落ち着いたあたりのクラスタ数が好ましい数だとみなす考えだそうです。
elbow <- function(dmatrix, kmax) { npts <- dim(dmatrix)[1] wss <- numeric(kmax) wss[1] <- (npts-1)*sum(apply(dmatrix, 2, var)) for(k in 2:kmax) { d <- dist(dmatrix, method="euclidean") pfit <- hclust(d, method="ward.D2") labels <- cutree(pfit, k=k) wss[k] <- wss.total(dmatrix, labels) } list(wss = wss) } elbow <- elbow(pmatrix,10) elbowframe <- data.frame(k=1:10,wss=elbow$wss) elbowframe <- melt(elbowframe, id.vars=c("k"), variable.name="measure", value.name="score") ggplot(elbowframe, aes(x=k, y=score, color=measure)) + geom_point(aes(shape=measure)) + geom_line(aes(linetype=measure)) + theme_dark() + scale_x_continuous(breaks=1:10, labels=1:10) + labs(x = "クラスタ数", y = "クラスタ内距離二乗和")
うーん。クラスタ数が5あたりからクラスタ内距離二乗和の値の減少幅は緩やかになってますが、この図からだと判断しずらくなんとも言えないです。
Calinski-Harabasz (CH) 基準（Pseudo F）
クラスタ内距離二乗和はクラスタ内のデータの凝集性に着目した指標であるのに対して、Calinski-Harabasz (CH) インデックス（Pseudo F）はクラスタ内の凝集性に加えてクラスタ間の離散性を考慮した指標となります。
# total sum of squares totss <- function(dmatrix) { grandmean <- apply(dmatrix, 2, FUN=mean) sum(apply(dmatrix, 1, FUN=function(row){sqr_edist(row, grandmean)})) } # wss Pseudo_F <- function(dmatrix, kmax) { npts <- dim(dmatrix)[1] # number of rows. totss <- totss(dmatrix) wss <- numeric(kmax) crit <- numeric(kmax) wss[1] <- (npts-1)*sum(apply(dmatrix, 2, var)) for(k in 2:kmax) { d <- dist(dmatrix, method="euclidean") pfit <- hclust(d, method="ward.D2") labels <- cutree(pfit, k=k) wss[k] <- wss.total(dmatrix, labels) } bss <- totss - wss # bssの計算 crit.num <- bss/(0:(kmax-1)) # bss/k-1 crit.denom <- wss/(npts - 1:kmax) # wss/n-k list(crit = crit.num/crit.denom, wss = wss, totss = totss, test= crit.num)# crit = (bss/k-1) / (wss/n-k) } # 視覚化 Fcrit <- Pseudo_F(pmatrix, 10) Fframe <- data.frame(k=1:10, ch=Fcrit$crit) Fframe <- melt(Fframe, id.vars=c("k"), variable.name="measure", value.name="score") ggplot(Fframe, aes(x=k, y=score, color=measure)) + geom_point(aes(shape=measure)) + geom_line(aes(linetype=measure)) + theme_dark() + scale_x_continuous(breaks=1:10, labels=1:10) + labs(x = "クラスタ数", y = "Pseudo F")
CH基準ではクラスタ数が2のときが好ましいという結果ができました。
実際にクラスタ数2でclusterboot関数を利用した結果は以下のとおりです。
cluster bootMeans 1 0.81 2 0.89
確かにJaccard係数の平均値（bootMeans）はともに0.8を超えており、すべてのクラスタの安定性はクラスタ数が5のときよりも高いといえるでしょう。しかし、デンドログラムを思い出すと分かるように、クラスタ数2では情報量が少なすぎて何も解釈できなくなってしまいます。
*1:詳細はテキストP.184を参照ください
