本次笔记内容：
cor.test()和cor()rcorr() {Hmisc}corr.test() {psych}
相关系数（correlation coefficient）用于描述两个变量之间的相关程度。一般在[-1, 1]之间。包括：
- pearson相关系数：适用于连续性变量，且变量服从正态分布的情况，为参数性的相关系数。
- spearman等相关系数：适用于连续性及分类型变量，为非参数性的相关系数。
在本次笔记中仅讨论连续型变量的相关系数。
# 示例数据有6个变量：
data("attitude")
head(attitude)
rating complaints privileges learning raises critical advance
1 43 51 30 39 61 92 45
2 63 64 51 54 63 73 47
3 71 70 68 69 76 86 48
4 61 63 45 47 54 84 35
5 81 78 56 66 71 83 47
6 43 55 49 44 54 49 34
cor.test()和cor()都是R自带包里的函数，两者差别仅为cor()只给出相关系数一个值，cor.test()给出相关系数，p值等。
你可以把数据的两组feature提出来进行相关性分析，看是否有相关性；也可以把包含多个feature的表格作为cor()input，得到的是一个对称的correlation matrix. 即所有feature两两比较的相关系数。然后你可以拿去各种可视化。cor.test()似乎不能这样用。
使用Hmisc包的rcorr()，可以得到correlation matrix的p值矩阵。当然rcorr()也可以像cor()那样，只计算两个feature之间的相关系数。
## 只把attitude中的rating和complaints作为input
cortest_ra_com <- cor.test(attitude$rating, attitude$complaints, method = "pearson")
cor_ra_com <- cor(attitude$rating, attitude$complaints, method = "pearson")
# 得到结果如下：
# > cortest_ra_com
# Pearson's product-moment
# correlation
# data: attitude$rating and attitude$complaints
# t = 7.737, df = 28,
# p-value = 1.988e-08
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
# 0.6620128 0.9139139
# sample estimates:
# cor
# 0.8254176
# > cor_ra_com
# [1] 0.8254176
## 把attitue中6个feature都作为input
cor_ <- cor(attitude, method = 'pearson')
View(cor_) # 如下图所示
library(Hmisc)
cortest <- rcorr(as.matrix(attitude), type = "pearson")
View(cortest$P) # 如下图所示
View(cor_)与View(cortest$r)结果一样
View(cortest$P)
如果你想比较attitude6个feature中前3个与后3个的关联，并且需要进行多重矫正，需要使用psych包的corr.test()。
你有关于一套sample的两套feature，比方说两个dataframe, 其行是相同的(sample)，列为不同的feature.那么可以corr.test(df1, df2, method= ...)来计算两组feature的相关系数并加以矫正。这时得到的output不是对称的，而是ncol(df1) * ncol(df2)
需要注意如果input为两个dataframe, 两者的row必须长度和顺序都一致。
library(psych)
cortest_psy <- corr.test(attitude[1:3], attitude[4:6], method = "pearson")
cortest_psy_sdj <- corr.test(attitude[1:3], attitude[4:6], method = "pearson", adjust = "fdr")
# 如果不矫正，即adjust ="none"，则其相关系数与P值其实和cor.test()等得到的一样。
可以根据P值，把P值做成*，**...这样的的significant levels，便于后面画热图。总的来说以下函数可以塞进去两个你想比较的dataframe，得到相关系数，矫正后的P值，校正后的P值significant levels矩阵，结合heatmap.2，就可以画图了...
R里做相关系数的函数茫茫多，不止这几个。以后如果要用到其他的再补上。
### fun_to_corr:
#### input: heat_in_1, heat_in_2: f_ra/g_ra/biochem/kegg_in...you can select the feature first
#### output: list(), t_cor is correlation index(-1~1), t_p is raw p-value, t_p_sig is formatted with significant levels
#### formatted significant levels: 0~0.001: **; 0.001~0.01: *; 0.01~0.1: +; >0.1 nothing
fun_to_corr <- function(heat_in_1, heat_in_2) {
t <- corr.test(heat_in_1, heat_in_2, use = "pairwise", method = "spearman", adjust = "fdr")
t_cor <- data.frame(t$r, check.names = FALSE)
t_p <- data.frame(t$p, check.names = FALSE)
cut_sig <- function(p) {
out <- cut(p, breaks = c(0, 0.001,0.01,0.1,1), include.lowest = T, labels = c("**", "*", "+", ""))
return(out)
}
t_p_sig <- apply(t_p, 2, cut_sig)
rownames(t_p_sig) <- rownames(t_p)
return(list(t_cor = t_cor, t_p_sig = t_p_sig, t_p = t_p))
}
Ref:
更多见STHDA的教程
corr.test的文档
网友评论