- はじめに
- パーマー群島(南極大陸)のペンギンさんのデータを活用する
- hclust関数での階層的クラスタリング
- 相関関係に基づく距離計算と階層的クラスタリング
- 階層的クラスタリングの結果をcutree関数やrect.hclust関数でグルーピングする
- pvclustによるマルチスケールブートストラップと階層的クラスタリング
- heatmap.2による、階層的クラスタリングとヒートマップ
- まとめ
はじめに
階層的クラスタリングは、以前からまとめてみたかった話です。 グルグル廻って、ようやく、ここまで辿り着いた感じです。
この記事では、R言語での階層的クラスタリングの基本ということで、 主に、hclustとheatmap.2という、2つの関数を扱います。 けど、関数群をまとめていると、他のパッケージとかも色々と盛り沢山になってしまった。。。
まず、hclust関数
は、statsパッケージ内の関数の1つで、 『非類似度の集合に関する階層的クラスター分析とその分析方法』を提供する、 R業界で最も有名な、階層的クラスタリング関数の1つです。
他方、gplotsパッケージが提供する heatmap.2関数
は、
従来の階層的クラスター分析に加えて、 ヒートマップや密度データなどの付加情報も一緒に描ける便利な関数となっています。
今回、ペンギンさんのデータを対象に、階層的クラスタリングをやってみることにします。
パーマー群島(南極大陸)のペンギンさんのデータを活用する
ペンギンさんのデータは、
palmerpenguinsパッケージ
で提供されている、 パーマーペンギンズ(palmerpenguins)のデータのことである*1。
このデータセットでは、 「南極のパーマー基地周辺のパーマー群島の島々で観察された アデリー、ヒゲペンギン、ジェンツーペンギンの成鳥のサイズ測定、クラッチ観察、 および血液同位体比のデータが掲載されている。 このデータは、Dr. Kristen Gorman氏とthe Palmer Station Long Term Ecological Research (LTER) Programによって収集され利用可能」となっています。
このデータそのままだと欠損値があるので、 まずは、それらを除去する作業とかをします。
#ペンギンさんのデータのロード install.packages("palmerpenguins") library(palmerpenguins) #ヘッダー表示 head(penguins) # A tibble: 6 × 8 # species island bill_length_mm bill_depth_mm # <fct> <fct> <dbl> <dbl> #1 Adelie Torgersen 39.1 18.7 #2 Adelie Torgersen 39.5 17.4 #3 Adelie Torgersen 40.3 18 #4 Adelie Torgersen NA NA #5 Adelie Torgersen 36.7 19.3 #6 Adelie Torgersen 39.3 20.6 # … with 4 more variables: flipper_length_mm <int>, # body_mass_g <int>, sex <fct>, year <int> #データフレームに変換 penguinsDF <- data.frame(penguins) #欠損数のカウント sum(is.na(penguinsDF)) #[1] 19 #欠損を除く penguinsDF01 <- na.omit(penguinsDF[,c("species", "island", "bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g")]) #欠損数のカウント sum(is.na(penguinsDF01)) #[1] 0 #speciesisland列の追加: species + island + 通し番号 penguinsDF01$species.island <- paste(penguinsDF01$species, penguinsDF01$island, 1:nrow(penguinsDF01)) penguinsDF02 <- penguinsDF01[,-c(1:2)] #最終的なデータセット head(penguinsDF02) # bill_length_mm bill_depth_mm flipper_length_mm body_mass_g species.island #1 39.1 18.7 181 3750 Adelie Torgersen 1 #2 39.5 17.4 186 3800 Adelie Torgersen 2 #3 40.3 18.0 195 3250 Adelie Torgersen 3 #5 36.7 19.3 193 3450 Adelie Torgersen 4 #6 39.3 20.6 190 3650 Adelie Torgersen 5 #7 38.9 17.8 181 3625 Adelie Torgersen 6
データとしては、欠損値(NA)を除いて、 "bill_length_mm"、"bill_depth_mm"、 "flipper_length_mm"、"body_mass_g"カラムのデータを使用しています。
hclust関数での階層的クラスタリング
「penguinsDF02」のデータフレームを加工して、数値データのみにします。 次に、scale関数で、ScalingとCentering補正する場合としない場合とで、 データフレームの作成を行なっています。
#データの加工 (数値データのみ) Column <- c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g") #補正無しの場合 penguinsDF03 <- penguinsDF02[,Column] #ScalingとCentering補正の場合 penguinsDF04 <- scale(penguinsDF03) head(penguinsDF04) # bill_length_mm bill_depth_mm flipper_length_mm body_mass_g #1 -0.8832047 0.7843001 -1.4162715 -0.5633167 #2 -0.8099390 0.1260033 -1.0606961 -0.5009690 #3 -0.6634077 0.4298326 -0.4206603 -1.1867934 #5 -1.3227986 1.0881294 -0.5628905 -0.9374027 #6 -0.8465718 1.7464261 -0.7762357 -0.6880121 #7 -0.9198375 0.3285561 -1.4162715 -0.7191859
これらのデータセットを使用して、 次に、hclustのagglomeration methodをいろいろと試してみました。 距離計算には、ユークリッド距離(Euclidean distance)を使っています。
ウォード法(ward.D)を使用する場合
補正有りのpenguinsDF04で、階層的クラスタリングをやってみました。
#"ward.D" method Xd <- dist(penguinsDF04, method="euclidean") h <- hclust(Xd, method = "ward.D") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
ちなみに、補正無しのpenguinsDF03でやって見ると。
#"ward.D" method Xd <- dist(penguinsDF03, method="euclidean") h <- hclust(Xd, method = "ward.D") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
この場合、補正無しなので、縦軸の高さがかなり大きくなります。
"ward.D" method + データフレームの転置
データの方向を変えてクラスタリングもできますが、 4変数を見るだけだと、あまり面白みはないかも。 t関数で、行列・データフレームの転置ができます。
#"ward.D" method + データフレームの転置 Xd <- dist(t(penguinsDF04), method="euclidean") h <- hclust(Xd, method = "ward.D") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(1, 5, 5, 1)) plot(h, hang=-1, cex=0.5, xlab = "", sub = "")
ウォードD2法(ward.D2)を使用する場合
#"ward.D2" method Xd <- dist(penguinsDF04, method="euclidean") h <- hclust(Xd, method = "ward.D2") #Or Xd^2も同じ結果 #h <- hclust(Xd^2, method="ward.D") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
私個人的には、デンドログラムのリーフの方の感じとしてはこれ形状が好きな感じです。
最近隣法(single)を使用する場合
#"single" method Xd <- dist(penguinsDF04, method="euclidean") h <- hclust(Xd, method = "single") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
群平均法(average)を使用する場合
#"average" method Xd <- dist(penguinsDF04, method="euclidean") h <- hclust(Xd, method = "average") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
McQuitty法(mcquitty)を使用する場合
#"mcquitty" method Xd <- dist(penguinsDF04, method="euclidean") h <- hclust(Xd, method = "mcquitty") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
メディアン法(median)を使用する場合
#"median" method Xd <- dist(penguinsDF04, method="euclidean") h <- hclust(Xd, method = "median") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
重心法(centroid)を使用する場合
#"centroid" method Xd <- dist(penguinsDF04, method="euclidean") h <- hclust(Xd, method = "centroid") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
私にとっては、"ward.D2" methodがデンドログラムに馴染みがあって、見易い気がします。
という事で、"ward.D2" methodに絞って、他の距離計算法についても試してみることにします。
相関関係に基づく距離計算と階層的クラスタリング
相関関係を使って距離計算を行い、階層的クラスタリングを行う場合を紹介します。
一般的に、相関係数は、「1」に近いほど相関が強い、つまりは『似ている』ことを示します。 他方、「0」に近いほど、「相関が弱い・相関関係が無い」ことを示します。
距離の計算上は、似ていない場合にお互いの距離が離れててほしいため、 『1 - 相関係数』を距離として定義することになります。 この距離関係では、『0』に近いほど、「似ている」ことを表せます。
スピアマン(Spearman)の順位相関係数と階層的クラスタリング
ここでは、『1 - (Spearmanの順位相関係数)』の関係を使って、 距離行列を作成して、階層的クラスタリングをやってみます。 また、これまでのdist()ではなく、代わりに、 as.dist()を使用するのがポイントとなります。
##1 - (Spearmanの順位相関係数)での距離行列 Xd <- as.dist(1 - cor(t(penguinsDF04), method = "spearman")) h <- hclust(Xd, method = "ward.D2") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
(補足として)
dist関数
は実際に距離行列を計算するのに対して、
as.dist関数
はオブジェクトを距離行列に強制的に変換するのみとのことです。
良く解らんけど・・・
https://lecture.ecc.u-tokyo.ac.jp/~aiwata/biostat_basic/2013/text4lec4_2.pdf

コサイン距離と階層的クラスタリング
コサイン距離(コサイン類似度)を使ったケースも紹介します。
コサイン類似度とは、内積空間の2つのベクトル間の類似度を表す尺度となります。
R内でのパッケージとしては、
LSA (Latent Semantic Analysis)パッケージのcosine()関数
を使います。
#lsaパッケージのインストール install.packages("lsa") library(lsa) ##1 - コサイン類似度での距離行列 Xd <- as.dist(1 - cosine(t(penguinsDF04))) h <- hclust(Xd, method = "ward.D2") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
ペンギンさんのデータでは、Spearman使うよりも、 コサイン距離の方が、階層クラスタリングとしてはキレイに見えますね。
混合データのときの距離行列計算と階層的クラスタリング
cluster::daisy()
を使用すると、
混合データ(数字以外の文字列のデータにも使える)でも
距離行列が計算できるようになります。
なんだか、便利そうな響きがします。
ここでは、文字列(Factorとなっているが)のデータが含まれている、 penguinsDF01のデータセットを使用します。
#clusterパッケージのインストール install.packages("cluster") library(cluster) #ディメンジョン表示 dim(penguinsDF01) #[1] 342 9 #詳細表示 str(penguinsDF01[,c(1:6)]) #'data.frame': 342 obs. of 6 variables: # $ species : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ... # $ island : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ... # $ bill_length_mm : num 39.1 39.5 40.3 36.7 39.3 38.9 39.2 34.1 42 37.8 ... # $ bill_depth_mm : num 18.7 17.4 18 19.3 20.6 17.8 19.6 18.1 20.2 17.1 ... # $ flipper_length_mm: int 181 186 195 193 190 181 195 193 190 186 ... # $ body_mass_g : int 3750 3800 3250 3450 3650 3625 4675 3475 4250 3300 ... #非類似度マトリクス計算 Xd <- cluster::daisy(penguinsDF01[,c(1:6)], metric="gower") str(Xd) #num [1:342, 1:342] 0 0.0447 0.0839 0.0742 0.069 ... # - attr(*, "dimnames")=List of 2 # ..$ : chr [1:342] "1" "2" "3" "5" ... # ..$ : chr [1:342] "1" "2" "3" "5" ... #クラスタリング h <- hclust(Xd, method = "ward.D2") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")

factoextra::get_dist()関数を使った相関計算と階層的クラスタリング
次に、factoextraパッケージのget_dist()関数
を使った方法を紹介します。
get_dist()関数
では、データ行列において行方向の距離行列を計算します。
先述した、cor関数使用の場合だとR上で計算式を組む必要がありますが、
factoextraパッケージのget_dist()関数
では、
標準でよく使われる、dist()関数では扱われていない、相関に基づく距離尺度をサポートしています。
そのため、methodとして、ピアソン(pearson)、スピアマン(spearman),ケンドール(kendall)の相関係数を選択できます。
#factoextraのインストール install.packages("factoextra") library(factoextra) ##ピアソン相関での距離行列 Xd <- factoextra::get_dist(x = penguinsDF04, method = "pearson") h <- hclust(Xd, method = "ward.D2") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "")
一応、『get_dist()とas.dist(1 - cor)の結果は同じか?』を検証してみました。
#pearson法で計算した場合 Xd1 <- factoextra::get_dist(x = penguinsDF04, method = "pearson") Xd2 <- as.dist(1 - cor(t(penguinsDF04), method = "pearson")) #差分計算 sum(Xd1 - Xd2) #[1] 0 #spearman法で計算した場合 Xd1 <- factoextra::get_dist(x = penguinsDF04, method = "spearman") Xd2 <- as.dist(1 - cor(t(penguinsDF04), method = "spearman")) #差分計算 sum(Xd1 - Xd2) #[1] 0
その結果、pearsonでもspearmanでも差分は「0」で同じ距離行列を返すことが分かりました。
階層的クラスタリングの結果をcutree関数やrect.hclust関数でグルーピングする
cutree関数やrect.hclust関数を使うと、 hclustで得られた階層的クラスタリングの樹木図を 希望するグループ数または高さ(距離)を指定して いくつかのグループに分割することができます。
#"ward.D2" method Xd <- dist(penguinsDF04, method="euclidean") h <- hclust(Xd, method = "ward.D2") #図示 par(family="HiraKakuProN-W3", mgp=c(2.5, 1, 0), xpd=T, mar=c(2, 5, 5, 1)) plot(h, hang=-1, cex=0.1, xlab = "", sub = "") #3グループに分割する hcut <- cutree(h, k = 3) #テーブル表示 table(hcut) #hcut # 1 2 3 #162 123 57 #クラスターを3つの四角で囲む rect.hclust(h , k = 3, border = 2:4) #高さhで切って、whichは番号(ツリーの左から右へ)でクラスタを選択する #デフォルトは which = 1:k です。 rect.hclust(h, h = 5, which = 7, border = "blue")
pvclustによるマルチスケールブートストラップと階層的クラスタリング
pvclustパッケージ
は、
階層型クラスター分析における不確実性を評価するためのパッケージです。
階層型クラスタリングの各クラスタについて、 マルチスケールブートストラップリサンプリングにより AU (Approximately Unbiased p-value) や BP (Bootstrap Probability)といった統計量をを算出できます。 クラスターのP値は、0 から 1 の間の値を示し、 そのクラスタがデータによってどの程度強く支持されるかを示します。
#pvclustパッケージのインストール install.packages("pvclust") library(pvclust) #マルチスケールブートストラップリサンプリング(non-parallel) Xd.pv <- pvclust(penguinsDF04, nboot=100, parallel=FALSE) #図示 plot(Xd.pv, cex.pv = 0.5, cex = 0.5, print.pv = c("au", "si", "bp"), hang=-1) pvrect(Xd.pv, pv="au")
図内で、
au
は、近似不バイアスp値(the approximately unbiased p-value)、
si
は、選択的推論のための近似不バイアスp値(apporximately unbiased p-values for selective inference)、
bp
は、ブーストラップ確率(boostrap probability)を示します。
t(penguinsDF04)
として、行と列を反転させて計算することもできて、
そうすると、変数が多くなってデンドログラムの解析っぽくなります。
heatmap.2による、階層的クラスタリングとヒートマップ
最後に、gplotsパッケージ内のheatmap.2関数を使用して、 拡張ヒートマップによる可視化をやってみます。 各種の引数は、以前に最適化したパラメータで行なっています。
#gplotsパッケージのインストール install.packages("gplots") library(gplots) #作図の設定 par(family="HiraKakuProN-W3", mgp=c(3, 1, 0), xpd=T) #5つのグループを事前に計算しておく Xd <- dist(penguinsDF04, method="euclidean") mat.c <- cutree(hclust(Xd, method = "ward.D2"), k=5) clust.col <- c("black", "red", "blue", "green", "yellow") #ヒートマップ実行 heatmap.2(penguinsDF04, col = colorpanel(100,"blue","white","red"), na.color="grey40", distfun = function(x) {dist(x, method="euclidean")}, hclustfun = function(x) {hclust(x, method="ward.D2")}, cexRow = 0.3, cexCol = 0.5, trace="none", symkey=FALSE, dendrogram = "both", Rowv = T, Colv = T, scale = "non", na.rm=T, RowSideColors=clust.col[mat.c], keysize = 0.7, densadj = 0.25, main="Enhanced Heat Map", key.par = list(cex=0.3, mar=c(4,4,3.5,0), mgp=c(2.5, 1, 0)), sepcolor="grey40",sepwidth=c(0.1,0.1), offsetRow = 0.5, offsetCol = 0.25, adjRow = c(0,0.25), adjCol = c(NA,0.25))

まとめ
階層的クラスタリングについて、色々とまとめてみました。
データ解析のチュートリアルとして使える機会があれば幸いです。
また、次回のクラスタリングの企画では、 dendextendやfactoextraパッケージを使った事例を紹介したいと思っています。
関連資料
https://www.ic.nanzan-u.ac.jp/~urakami/pdf/19Rpdf/19m_23.pdf
*1:先日のTokyo.Rで見かけたので、このデータを使いたくなった