Rで、サンプルデータ作製からヒートマップ表示とクラスタリング表示 (Afiinity propagation含む) まで

Rで、サンプルデータ作製して、ヒートマップ表示と、k-medoids (pam) と Afiinity propagation含むいくつかのクラスタリング表示を作る機会がありました。せっかくなので、記載します。今回ようなデータの処理と目的では、k-means のような方法が適していました。DBSCAN などは他をご参考をいただけますようお願いします。ひたすら駆け足で、忘備録程度の内容です。

対象は、50人分、10種類の観測値のサンプルデータです。

目的は、ヒトごと、10種類の観測値ごとのそれぞれのクラスターを見付けることです。

平均値とSD の異なる正規分布で、50人分、10種類の観測値のサンプルデータを作ります。

A1=cbind(
  matrix(rnorm(50, -1, 1),nc=5), 
  matrix(rnorm(50,0,1),nc=5)
  )
A2=cbind(
  matrix(rnorm(30,-0.5,1),nc=3), 
  matrix(rnorm(50,0,2),nc=5),
  matrix(rnorm(20,-3,1),nc=2)
  )
A3=cbind(
  matrix(rnorm(20,0.5,2),nc=2),
  matrix(rnorm(30,-2,1),nc=3),
  matrix(rnorm(20,0.2,0.2),nc=2),
  matrix(rnorm(30,1,1),nc=3)
  )

A4=cbind(
  matrix(rnorm(10,2,0.8),nc=1),
  matrix(rnorm(40,1,0.7),nc=4),
  matrix(rnorm(10,-0.5,0.1),nc=1),
  matrix(rnorm(40,0,1),nc=4)
  )
A5=cbind(
  matrix(rnorm(10,-0.1,0.2),nc=1),
  matrix(rnorm(40,-2,1),nc=4),
  matrix(rnorm(10,2,1),nc=1),
  matrix(rnorm(40,0.4,1),nc=4)
  )
M_Example = rbind(A1,A2,A3,A4,A5)

データを50人の個人をランダムに並べ変えます。 クラスタライングの為には、正規化しておく方が良いことが多いです。 10種の観測値ごとに正規化し、data に納めます。 ここまでを1行で、次からは、行列に名前を付けます。

data = scale((M_Example[order(runif(50)),]))

col_data =c()
for (i in 1:10){
  col_data <- c(col_data, paste("data", as.character(i), sep=""))
}

row_data =c(  )
for (i in 1:50){
  row_data <- c(row_data, paste("person", as.character(i), sep=""))
}

colnames(data)<-col_data
rownames(data)<-row_data

できあがったサンプルデータを、順番を変えずに、ヒートマップで示します。ヒートマップは、一貫して heatmap.2 を使いました。

library(gplots)

# 順番を変えずに描いてみましょう。
Colour <- colorRampPalette( c("blue3","white","red3") )
heatmap.2(data,
  dendrogram = "none",
  Colv=FALSE,
  Rowv=FALSE,
  scale="none",
  density.info="none",
  col=Colour(64),
  main="Heatmap (no cluster)", margin=c(8,10),
  trace = "none"
)

平均値を白、低いものを青、高いものを赤で示しています。以下ではスケールバーは、同じなので、省略します。

つづいて、デフォルトであるところの最長距離法 (furthest neighbor method) または 完全連結法 (complete linkage method) と呼ばれる方法によるもの。

Colour <- colorRampPalette( c("blue3","white","red3") )
heatmap.2(data,
  dendrogram = "both",
  Colv=T,
  Rowv=T,
  scale="none",
  density.info="none",
  col=Colour(256),
  main="Heatmap (complet)",
  trace = "none"
)

つついて、Ward 法を書きます。Dist と hclust でクラスタリングを別にする方が、中身が見えるので、やりやすいでしょう。

library(amap)
# 距離関数で類似度のアルゴリズムとしてCosine係数を使用する
# (method="pearson")は(1 - cosine)と同義である
d1<-Dist(data, method="pearson")
d2<-Dist(t(data), method="pearson")

# クラスタリング
c1<-hclust(d1, method="ward.D")
c2<-hclust(d2, method="ward.D")

Colour <- colorRampPalette( c("blue3","white","red3") )
heatmap.2(data,
  dendrogram = "both",
  Colv=as.dendrogram(c2),
  Rowv=as.dendrogram(c1),
  scale="none",
  key = T,
  col=Colour,
  density.info="none",
  main="Heatmap (Ward algorythm)",
  trace = "none"
)

k-medoids です。パッケージユーザーの水準なら、k-medoids は k-means の改良版と考えて問題が生じることは少いのではないでしょうか。k-means は、走らせる度に、クラスターが変化しますが、k-medoids は変化しません。これだけでも、パッケージユーザーには利点があると考えます。今回は独断と偏見で、両方向ともクラスタラリング数は、8 としています。樹形図は出せないので、色分けでクラスタラリングを表現しています。

# クラスタリング
library(cluster)
# mat_k4_th_norm <- t(scale(t(mat_k4_th)))
cluster_n  <- 8
cluster_nc <- 6
result_pam <- pam(data, cluster_n)
pam_clust <- result_pam$cluster

result_pam_c <- pam(t(data), cluster_nc)
pam_clust_c <- result_pam_c$cluster

df_heat <- data[order(pam_clust, data[,1]),order(pam_clust_c)]


# Heatmap の出力; 
library(gplots)
library('RColorBrewer')


# 右横にクラスター帯を作るための色を決めます。
col_col=c()
i=1
while ( i < cluster_n    ){
col_col=c(col_col,
  rep(brewer.pal(8,"Accent")[1],t(table(pam_clust))[i]),
  rep(brewer.pal(8,"Accent")[2],t(table(pam_clust))[i+1]),
  rep(brewer.pal(8,"Accent")[3],t(table(pam_clust))[i+2]),
  rep(brewer.pal(8,"Accent")[4],t(table(pam_clust))[i+3]),
  rep(brewer.pal(8,"Accent")[5],t(table(pam_clust))[i+4]),
  rep(brewer.pal(8,"Accent")[6],t(table(pam_clust))[i+5]),
  rep(brewer.pal(8,"Accent")[7],t(table(pam_clust))[i+6]),
  rep(brewer.pal(8,"Accent")[8],t(table(pam_clust))[i+7])
  )
  i= i + 6
}

#上の帯
col_col_c=c()
i=1
while ( i < cluster_nc    ){
col_col_c=c(col_col_c,
  rep(brewer.pal(8,"Accent")[1],t(table(pam_clust_c))[i]),
  rep(brewer.pal(8,"Accent")[2],t(table(pam_clust_c))[i+1]),
  rep(brewer.pal(8,"Accent")[3],t(table(pam_clust_c))[i+2]),
  rep(brewer.pal(8,"Accent")[4],t(table(pam_clust_c))[i+3]),
  rep(brewer.pal(8,"Accent")[5],t(table(pam_clust_c))[i+4]),
  rep(brewer.pal(8,"Accent")[6],t(table(pam_clust_c))[i+5])
  )
  i= i + 6
}


# ヒートマップ自身の色
Colour <- colorRampPalette( c("blue3","white","red3") )

heatmap.2(as.matrix(df_heat),
  dendrogram = "none",
  Rowv = "none",
  Colv = "none",
  scale="none",
  col= Colour(64),
  main="Heatmap (pam)",
  trace = "none",
  RowSideColors= (t(col_col[1:dim(df_heat)[1]])),
  ColSideColors= (t(col_col_c[1:dim(df_heat)[2]]))
)

最後に、affinity propagation です。個人的には、affinity propagation が好きです。帯でクラスタリングを表現しています。この帯を描くところは、一般化が不十分のあコードです。クラスタリングの数によって、書き変えないと走らないことがあるので、注意して下さい。"q =" のオプション (0-1) で、クラスタリングの数をある程度コントロールできます。たぶん、樹形図の高さに相当していると考えます。

library(apcluster)

result_apc <- apcluster(negDistMat(r=2), data)
label_apc <- labels(result_apc,type="enum")
cluster_n <- length(table(label_apc))


result_apc_c <- apcluster(negDistMat(r=2), t(data), q=0.1)
label_apc_c <- labels(result_apc_c,type="enum")
cluster_nc <- length(table(label_apc_c))

df_heat <- data[order(label_apc, data[,1]),order(label_apc_c)]

# Heatmap の出力; 
library(gplots)
library('RColorBrewer')


# 右横にクラスター帯を作るための色を決めます。
col_col=c()
i=1
while ( i < cluster_n    ){
col_col=c(col_col,
  rep(brewer.pal(8,"Accent")[1],t(table(label_apc))[i]),
  rep(brewer.pal(8,"Accent")[2],t(table(label_apc))[i+1]),
  rep(brewer.pal(8,"Accent")[3],t(table(label_apc))[i+2]),
  rep(brewer.pal(8,"Accent")[4],t(table(label_apc))[i+3]),
  rep(brewer.pal(8,"Accent")[5],t(table(label_apc))[i+4]),
  rep(brewer.pal(8,"Accent")[6],t(table(label_apc))[i+5]),
  rep(brewer.pal(8,"Accent")[7],t(table(label_apc))[i+6]),
  rep(brewer.pal(8,"Accent")[8],t(table(label_apc))[i+7])
  )
  i= i + 6
}

col_col_c=c()
i=1
while ( i < cluster_nc    ){
col_col_c=c(col_col_c,
  rep(brewer.pal(8,"Accent")[1],t(table(label_apc_c))[i]),
  rep(brewer.pal(8,"Accent")[2],t(table(label_apc_c))[i+1])
  )
  i= i + 3
}


# ヒートマップ自身の色
Colour <- colorRampPalette( c("blue3","white","red3") )

heatmap.2(as.matrix(df_heat),
  dendrogram = "none",
  Rowv = "none",
  Colv = "none",
  scale="none",
  col= Colour(64),
  main="Heatmap (Affinity propagation)",
  trace = "none",
  RowSideColors= (t(col_col[1:dim(df_heat)[1]])),
  ColSideColors= (t(col_col_c[1:dim(df_heat)[2]]))
)

クラスタリングの内訳は、heatmap() の中に結果を入れると描いてくれます。

library(viridis)
heatmap(result_apc, col=magma(64))
heatmap(result_apc_c, col=magma(64))

スケールは、省略で恐縮ですが、viridis パッケージの magma を使っています。明確な平均を0としていないような非正規化データの場合は、magma などの色の方が、差が見えます。

以上、ひたすら駆け足ですが、おわりです。