原文链接
相关性分析的应用场景
一些样本,每个样本会测一些指标,我想初步探索一下这些指标之间是否存在关联。具体场景:我收集了好多个品种的苹果成熟果实,每个品种的苹果我都会测一些指标,比如表型指标:果重;生理指标:可溶性糖,有机酸,花青素含量等等。
做完实验数据整理到excel中,另存为csv格式
image.png
数据是我胡编乱造的,没有实际意义!
读入数据
csvpath<-file.choose()
csvpath
df<-read.csv(csvpath,header=T,row.names = 1)
df
这样就把数据读进来存储到df里了
R语言里自带的相关性分析的函数是cor()
,直接将数据放到括号里就可以了。默认的皮尔逊相关性分析
> cor(df)
fruit_weight soluble_sugar organic_acid anthocyanin
fruit_weight 1.00000000 0.06342157 -0.2647533 0.1038605
soluble_sugar 0.06342157 1.00000000 0.2580373 -0.2590438
organic_acid -0.26475334 0.25803726 1.0000000 -0.2241183
anthocyanin 0.10386047 -0.25904381 -0.2241183 1.0000000
通过method参数指定其他方法
> cor(df,method = 'sperman')
Error in match.arg(method) :
'arg' should be one of “pearson”, “kendall”, “spearman”
> cor(df,method = 'spearman')
fruit_weight soluble_sugar organic_acid anthocyanin
fruit_weight 1.0000000 0.1357143 -0.1714286 0.1892857
soluble_sugar 0.1357143 1.0000000 0.2821429 -0.2000000
organic_acid -0.1714286 0.2821429 1.0000000 -0.2142857
anthocyanin 0.1892857 -0.2000000 -0.2142857 1.0000000
但是论文里的相关性分析通常都是带有p值,就是右上角会有星号。可以借助Hmisc
包中的rcorr
函数
这个函数要求的输入数据格式是矩阵,同过csv文件读入的数据格式是数据框,需要借助函数as.matrix()
进行转换
library(Hmisc)
res2<-rcorr(as.matrix(df))
运行完以后res2里面存储3个内容,可以通过$符号获取
> res2$r
fruit_weight soluble_sugar organic_acid anthocyanin
fruit_weight 1.00000000 0.06342157 -0.2647533 0.1038605
soluble_sugar 0.06342157 1.00000000 0.2580373 -0.2590438
organic_acid -0.26475334 0.25803726 1.0000000 -0.2241183
anthocyanin 0.10386047 -0.25904381 -0.2241183 1.0000000
> res2$n
fruit_weight soluble_sugar organic_acid anthocyanin
fruit_weight 15 15 15 15
soluble_sugar 15 15 15 15
organic_acid 15 15 15 15
anthocyanin 15 15 15 15
> res2$P
fruit_weight soluble_sugar organic_acid anthocyanin
fruit_weight NA 0.8223325 0.3402882 0.7126110
soluble_sugar 0.8223325 NA 0.3531301 0.3511885
organic_acid 0.3402882 0.3531301 NA 0.4219767
anthocyanin 0.7126110 0.3511885 0.4219767 NA
r是相关性系数,n是样本个数,p是相关性检验的p值
接下来我想看看谁跟谁的相关性比较高,比如筛选相关系数绝对值大于0.8。矩阵筛选我还不知道如何实现。原文自己写了一个函数,将矩阵转换为数据框,这样筛选起来就容易很多了。
函数是
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]
)
}
两个参数,一个是相关性,一个是p值
> flattenCorrMatrix(res2$r,res2$P)
row column cor p
1 fruit_weight soluble_sugar 0.06342157 0.8223325
2 fruit_weight organic_acid -0.26475334 0.3402882
3 soluble_sugar organic_acid 0.25803726 0.3531301
4 fruit_weight anthocyanin 0.10386047 0.7126110
5 soluble_sugar anthocyanin -0.25904381 0.3511885
6 organic_acid anthocyanin -0.22411828 0.4219767
筛选一个相关系数绝对值大于0.25的
> df1<-flattenCorrMatrix(res2$r,res2$P)
> abs(df1$cor)>0.25
[1] FALSE TRUE TRUE FALSE TRUE FALSE
> df1[abs(df1$cor)>0.25,]
row column cor p
2 fruit_weight organic_acid -0.2647533 0.3402882
3 soluble_sugar organic_acid 0.2580373 0.3531301
5 soluble_sugar anthocyanin -0.2590438 0.3511885
接下来就是数据展示了,一个是表格,一个是图。
接下来介绍画图:
一种展示方法
library(corrplot)
corrplot(res2$r,type="upper",tl.col ="black",tl.srt = 45)
image.png
另外一种展示方法
install.packages("PerformanceAnalytics")
library(PerformanceAnalytics)
chart.Correlation(df,histogram = T,pch=19)
image.png
还可以选择用热图来展示
col<-colorRampPalette(c("blue","white","red"))(20)
heatmap(x=res2$r,col=col,symm=T)
heatmap(x=res2$r,col=col,symm=F)
image.png
暂时还不知道symm参数的作用是啥?
数据大家完全可以自己构造,原文用到的数据是R本身自带的例子mtcars,但是各项指标可能不太好理解。所以我就自己随便伪造了一份数据。如果大家需要可以给我公众号留言
欢迎大家关注我的公众号
小明的数据分析笔记本
网友评论