- はじめに
- パーマー群島（南極大陸）のペンギンさんのデータを活用する
- 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で見かけたので、このデータを使いたくなった