Rのgraphics::heatmapのクラスター分析（樹形図）には罠がある

Rのgraphics::heatmapには予想し難い挙動があります。ヘルプや解説ページでも言及されていないようなので指摘しておきたいと思います。

観測数を集計した、行で標本、列で属性を表す行列を用意します。

M <- matrix(c(1,6,4,2,5,8,1,6,4), 3)
rownames(M) <- sprintf("sample #%d", 1:3)
colnames(M) <- sprintf("factor %s", letters[1:3])

頻度の集計なのでヒートマップがプロットできます。やってみましょう。

heatmap(M, Colv=NA, margin = c(5, 10))

heatmap（デフォルト）

ヒートマップと樹形図（デンドログラム）の対応が変です。sample #1と#2より、sample #1と#3の方が近いはずですが、そうは分析されていません。

これはヒートマップとクラスター分析で使う数字が異なるためです。graphics::heatmapはデフォルトで行ごとに標準化した数値でヒートマップの描画をします。属性ごとの観測数の比率が同じsample #1と#3は同じ配色になります。一方、graphics::heatmapはデフォルトで生の数字でクラスター分析をかけて、樹形図を作成します。

解決方法としては、scale = "none"をつけて、ヒートマップの表示をクラスター分析にあわせるか、以下のようにクラスター分析をカスタマイズする必要があります。

heatmap(M, Rowv = NULL, Colv = NA, 
  dist = \(X){
    margin <- 1 # 1だと行ごと，2だと列ごと
    m <- apply(X, margin, mean) # 行ごとの平均
    s <- apply(X, margin, sd) # 行ごとの標準偏差
    sX <- apply(X, c(2, 1)[margin], \(x) (x-m)/s) # 標準化
    sX[is.na(sX)] <- 0 # 標準偏差0でNAのときは0に置換
    if(2 == margin) sX = t(sX)
    dist(sX)
  }, hclustfun = \(d){
    hclust(d, method = "complete")
  },
  scale = "row", # デフォルト
  margins = c(5, 10))

heatmap（パラメーターを細かく指定）

たまたま説明*1のために3×3で挙動を確認したので気づいたのですが、これに気づかず使っている人がいそうでちょっと怖いですね。

*1:Rで頻度をプロット（hist/heatmap/mosaic/cd）