1. 矩陣相關性計算方法
base::cor/cor.test
R基礎函數cor或cor.test都可計算相關性系數，但cor可直接計算矩陣的相關性，而cor.test不可。
兩者計算非矩陣時，cor僅得到相關系數，而cor.test還能得到pvalue。
library(ggplot2)
cor(mtcars)
cor.test(mtcars) #error
cor.test(mtcars,mtcars) #error
cor(mtcars$mpg,mtcars$cyl) #only cor
x=cor.test(mtcars$mpg,mtcars$cyl) #cor and pvalue
x$estimate
x$p.value
可以用基礎函數cor得到相關性矩陣，再自己編寫腳本獲得pvalue矩陣。
M = cor(mtcars)
#自編寫函數得到pvalue矩陣
cor.mtest <- function(mat, ...) {
mat <- as.matrix(mat)
n <- ncol(mat)
p.mat<- matrix(NA, n, n)
diag(p.mat) <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
tmp <- cor.test(mat[, i], mat[, j], ...)
p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
}
}
colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
p.mat
}
matrix_p=cor.mtest(mtcars)
psych::corr.test
使用psych包中的corr.test函數，可直接獲得矩陣相關性系數和pvalue（也可用於非矩陣），而且還可直接得到矯正后的pvalue。
library(psych)
corr.test(mtcars)
cor <- corr.test(mtcars,
method = "pearson",
adjust = "fdr") #同p.adjust函數
cor$r
cor$p
cor$p.adj #但得到的是向量，數目也不對
test <- p.adjust(cor$p,method = "fdr")
identical(cor$p.adj,test) #不等
Hmisc::rcorr
使用Hmisc包中的rcorr函數，直接得到相關性系數和pvalue矩陣。
library(Hmisc)
#注意要將數據框轉換為矩陣
cor.mat <- rcorr(as.matrix(mtcars), type = "pearson")
cor.mat$r
cor.mat$P
可視化時，pvalue矩陣對角線的顯著性我們不必要展示，可以替換下。另外，如果后續不展示全部矩陣，只展示過了設置條件的部分，則可進行過濾。
# # only keep comparisons that have some abs. correlation >= .5 (optional)
# keep <- rownames(cor.mat$r)[rowSums(abs(cor.mat$r)>=0.5) > 1]
# cor.mat <- lapply(cor.mat, function(x) x[keep, keep])
# set diagonal to 1, since it is not interesting and should not be marked
diag(cor.mat$P) <- 1
其他工具
其他還有工具，如ggcor + ggcorrplot， 但不建議使用，增加學習成本，以上方法足以成對所有情況。
另外統計和繪圖R包rstatix也可計算相關矩陣，顯示和標記顯著性水平，而且可以gather和spread相關性矩陣，可tidyverse語法類似。這個包值得好好學習：https://rpkgs.datanovia.com/rstatix/index.html
2. 相關性矩陣轉化為兩兩相關
一般來說，我們得到的是相關性系數矩陣和pvalue矩陣，但輸出數據時最好轉換為兩兩之間的行列式格式。
這種轉換以上的rstatix包可輕松解決。
請參考：https://rpkgs.datanovia.com/rstatix/reference/cor_reshape.html
另外，我們也可自己編寫腳本得到：
flattenCorrMatrix <- function(cormat, pmat) {
ut <- upper.tri(cormat)
data.frame(
row = rownames(cormat)[row(cormat)[ut]],
column = rownames(cormat)[col(cormat)[ut]],
cor =(cormat)[ut],
p = pmat[ut]
)
}
res <- flattenCorrMatrix(cor.mat$r, cor.mat$P)
res
3. 可視化
得到了相關性和pvalue兩個矩陣，我們一般以熱圖展示為好。
corrplot
經典的相關性展示工具。很多可選樣式：https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html
我僅展示幾個案例，更多參數自己調節。
#僅cor
corrplot.mixed(M)
#cor,僅0.05
corrplot.mixed(M,
insig = 'label_sig',
p.mat=matrix_p,
pch.cex = 0.9,
pch.col = 'grey20')
#細分
corrplot(M,
p.mat = matrix_p,
tl.pos = 'd',
order = 'hclust',
type = "upper",
#addrect = 2,
insig = 'label_sig',
sig.level = c(0.001, 0.01, 0.05),
pch.cex = 0.9,
pch.col = 'grey20')
gplots::heatmap.2
相對於上圖，我更喜歡用熱圖來展示。
library(RColorBrewer)
library(gplots)
my_palette <- colorRampPalette(c("blue","white","red")) (100)
# plot heatmap and mark cells with abs(r) >= .5 and p < 0.05
heatmap.2(cor.mat$r,
# cexRow = .35, cexCol = .35,
trace = 'none',
# key.title = 'Spearman correlation',
# keysize = .5, key.par = list(cex=.4),
notecol = 'black', srtCol = 30,
col = my_palette,
cellnote = ifelse(cor.mat$P < 0.05 & abs(cor.mat$r)>=0.5, "*", ""))
以上我僅標出相關性絕對值大於0.5，pvalue<0.05的數據。當然可以做更細致划分。
pheatmap
pheatmap參數更好調些，看個人喜好。
#pheatmap
pheatmap(cor.mat$r,
color = my_palette,
display_numbers = ifelse(cor.mat$P < 0.05 & abs(cor.mat$r)>=0.5, "*", ""))
Ref:
https://www.jianshu.com/p/b76f09aacd9c
https://chowdera.com/2020/12/20201218185101270B.html
https://stackoverflow.com/questions/66305232/r-how-to-plot-a-heatmap-that-shows-significant-correlations
http://www.sthda.com/english/wiki/correlation-matrix-an-r-function-to-do-all-you-need
http://www.sthda.com/english/wiki/correlation-matrix-a-quick-start-guide-to-analyze-format-and-visualize-a-correlation-matrix-using-r-software