统计与科研生物信息学

R语言相关性分析简单小例子

2020-04-10  本文已影响0人  小明的数据分析笔记本

原文链接

http://www.sthda.com/english/wiki/correlation-matrix-a-quick-start-guide-to-analyze-format-and-visualize-a-correlation-matrix-using-r-software#at_pco=smlwn-1.0&at_si=5e8f19ae4cd478e7&at_ab=per-2&at_pos=0&at_tot=1

相关性分析的应用场景

一些样本,每个样本会测一些指标,我想初步探索一下这些指标之间是否存在关联。具体场景:我收集了好多个品种的苹果成熟果实,每个品种的苹果我都会测一些指标,比如表型指标:果重;生理指标:可溶性糖,有机酸,花青素含量等等。

做完实验数据整理到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,但是各项指标可能不太好理解。所以我就自己随便伪造了一份数据。如果大家需要可以给我公众号留言

欢迎大家关注我的公众号
小明的数据分析笔记本

公众号二维码.jpg
上一篇下一篇

猜你喜欢

热点阅读