[R][pheatmap] 有意差アスタリスクをヒートマップ内に表示

ライブラリ

library(pheatmap)
library(psych)

各変数間の相関係数をヒートマップで表す

相関係数を算出するcor関数をdata.frameに対して行うと、各列間の相関係数が算出できます。
mtcarsデータを使って実行してみます。

cor(mtcars)

この相関係数matrixをpheatmapで図示すると相関関係をわかりやすく表現できます。

pheatmap(cor(mtcars))

相関検定

相関検定を行うcor.test関数は2つのベクトルにしか対応していません。

data.frameの各列間の相関検定を行うにはpsychパッケージのcorr.test関数が使用できます。
もしくはppcorパッケージのpcor関数で偏相関係数の検定をしても良いでしょう。(ここでは解説しません。)

corr.testの引数確認

args(corr.test)

デフォルトではpearson相関検定をpairwiseに行い、多重検定補正はholm法で行います。

corr.testの実行

hoge <- corr.test(mtcars)

corr.testの返り値はリストです。
rに相関係数、p.adjに調整済みp値が保存されています。

corr.test出力の中身

pheatmap用にp値のmatrixを成型

この調整済みp値を使ってヒートマップ上に有意水準の*を入れていきたいのですが、上三角行列の値のベクトルになっており、調整済みp値は (行数 * 列数 - 対角成分数)/2 の数しかありません。

まずは調整済みp値を本来の行列サイズにしていきます。

調整済みp値のmatrix作成

# 同じサイズのNAのmatrixを用意
tmp <- matrix(data = NA, nrow = nrow(hoge$p), ncol = ncol(hoge$p))
# 上三角行列に調整済みp値を代入
tmp[upper.tri(tmp)] <- hoge$p.adj

# 上三角行列を転置して、下三角行列に変換 (上三角行列の箇所はNA)
p <- t(tmp)
# 上三角行列のNAの箇所に調整済みp値を代入
p[upper.tri(p)] <- hoge$p.adj

完成した調整済みp値のmatrix
ちなみに、corr.testの返り値pの行列は上三角が調整済みp値で下三角が調整前のp値なので、以下のようにしても同様の結果になります。

p値matrixの下三角に上三角の値を入れる

p <- t(hoge$p)
p[upper.tri(p)] <- hoge$p[upper.tri(hoge$p)]

次にp値の値に応じて*や空白セルに変換していきます。

値を記号に変換

for(i in 1:nrow(p)){
    for(j in 1:ncol(p)){
        tmp <- p[i,j]
        if(is.na(tmp)){
            symbol <- ""
        } else if(tmp < 0.001){
            symbol <- "***"
        } else if(tmp < 0.01){
            symbol <- "**"
        } else if(tmp < 0.05){
            symbol <- "*"
        } else{
            symbol <- ""
        }
        p[i,j] <- symbol
    }
}

最後にpheatmapのdisplay_numbers引数に上記で作成したmatrixを指定して完成です。

display_numbers引数

pheatmap(hoge$r, display_numbers = p)