Rによる階層的クラスタ分析の備忘録

教師なし学習の代表的な手法の1つであるクラスタ分析について個人的な論点整理を兼ねてZumel氏とMount氏共著の"Practical Data Science with R"の第8章をまとめてみます(端折ったり本文にないコードを追加したりしてます)。
ちなみに、第8章は以下の出版社サイトからFREEで入手可能です。

Manning | Practical Data Science with R


目次

データ

  • 欧州25ヵ国における9つの食品群(赤身肉、白身肉、卵、牛乳、魚、穀物、でんぷん質食品、豆類・堅果類・油種、果物・野菜)のタンパク質摂取量を測定したデータであり、著者のGitHubから入手可能です(Proteinというファイル名)

github.com

  • このようなデータセットを使用し、クラスタ分析により、食品ごとのたんぱく質摂取量のパターンにもとづいて25ヵ国をグルーピングしてみます。

使用するパッケージ

> 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つのクラスタに分類

f:id:joure:20170427222125p:plain
上図より、赤い枠で囲まれた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)

f:id:joure:20170429001315p:plain
クラスタ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より大きな値を有するクラスタはとても安定している(意味のあるクラスタを形成している)と判断されるらしい。上記の結果をみると、

ということが分かります。
デンドログラムからクラスタ数を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 = "クラスタ内距離二乗和")

f:id:joure:20170429162906p:plain

うーん。クラスタ数が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") 

f:id:joure:20170429221237p:plain

CH基準ではクラスタ数が2のときが好ましいという結果ができました。
実際にクラスタ数2でclusterboot関数を利用した結果は以下のとおりです。

cluster	bootMeans
1	0.81
2	0.89

確かにJaccard係数の平均値(bootMeans)はともに0.8を超えており、すべてのクラスタの安定性はクラスタ数が5のときよりも高いといえるでしょう。しかし、デンドログラムを思い出すと分かるように、クラスタ数2では情報量が少なすぎて何も解釈できなくなってしまいます。

結論

クラスター数を決定する唯一の方法はなさそう。分析対象のドメイン知識などヒューリステックに決めることも大切だと思われる。

*1:詳細はテキストP.184を参照ください