|
相关性分析通常用来定量描述两个变量之间的联系，正相关？负相关？不存在相关？等。常见相关性指标：Pearson相关系数（参数检验），Spearman及Kendall（非参数检验）。本期与分享R语言相关分析时的小技巧：1) 一个矩阵中等长度两两变量的相关性，常用到cor() 函数和cor.test() 函数；2) 一个矩阵中非等长度两个变量的相关性，常用corr.test()函数；3) 两个矩阵中变量对应的相关性，这里编了一个函数与大家分享。 (本期不做偏相关分析的讲解）1）Pearson相关系数（积差相关系数）通常用来表示相关性大小中最常用的指标。cor系数介于-1~1之间，越接近0相关性越弱，越接近-1 or 1则相关性越强，并存在对应的负正相关。2）Spearman等级相关系数（秩相关系数）利用两个变量的秩次大小进行比较，属于非参数检验。3）Kendall等级相关系数（和谐系数）也属于一种非参数检验。等长度是两个变量的样本量相等。主要利用 cor() 函数 和 cor.test() 函数对不同方法的讲解。以一组满足正态分布的连续变量和不满足正态分布的有序变量进行对比。例子讲解前，我们先用？cor（）/ help（cor）了解一下函数的用法和命令。
函数的基本命令：cor(x, use= , method= )x为一个数值向量，矩阵或者数据框都可以； y这里默认为NULL； na.rm 表示缺失值的逻辑，计算比较时是否保留缺失值（T表示去除缺失值，但是这里只能用于var函数中）；
use 指定缺失数据的处理方式。可选项：all.obs(假设不存在缺失数据)、everything(数据存在缺失值时，相关系数计算结果会显示missing)、complete.obs(行删除)、pairwise.complete.obs(成对删除)； method 指定相关系数的类型。可选类型为pearson、spearman、kendall。
清楚了函数对应命令的参数后，进行例子的理解 利用str（）检查数据结构的结果，num表示数值型的连续变量。#1 构建一组符合正态分布的数据
set.seed(1)
a <- rnorm(n = 50, mean = 1, sd = 5)
b <- rnorm(50, 2, 3)
c <- rnorm(50, 3, 2)
d <- rnorm(50, 0, 1)
# 将上述构建的数据组合成数据框, 并检查数据结构（num数值型）
e <- data.frame(a, b, c, d);str(e)
# pearson
corr1 <- cor(e, method = "pearson");corr1
# spearman
corr2 <- cor(e, method = "spearman");corr2
# kendall
corr3 <- cor(e, method = "kendall");corr2输出结果：当使用参数与非参数的方法对符合正态分布的连续变量进行相关性比较时，结果是存在差异的。Pearson为一个结果，而Spearman和kendall为两一个结果。
Example 2 -- 不满足正态分布的有序变量 利用as.integer（）将向量转换成整数型的，使其变成有序变量，数据结构为int，表示整数型。 #2 构建一组不满足正态分布的数据
set.seed(2)
a1 <- as.integer(runif(n = 50, min = 1, max = 10))
b1 <- as.integer(runif(n = 50, min = 0, max = 3))
c1 <- as.integer(runif(n = 50, min = 1, max = 6))
d1 <- data.frame(a1, b1, c1); str(d1)
# pearson
corr.1 <- cor(d1, method = "pearson");corr.1
# spearman
corr.2 <- cor(d1, method = "spearman");corr.2
# kendall
corr.3 <- cor(d1, method = "kendall");corr.2输出结果：从结果也能发现不同方法进行的相关性分析确实存在差异，但是具体哪种更好，可能需要大家多做尝试。 
接下来，利用cor.test（）函数对上面例子做详细展示，可以看得到对应的p值大小，有利于大家解读和后续可视化。还是先咨询下该函数的用法和对应参数。 
函数的基本命令：cor.test(x, y, alternative= , method= ,...) 
这里不做详细介绍，有问题的话可以评论留言~但这里需要注意的是：不能直接用数据框进行分析，而是需要用两两向量，因此用$提取对应向量进行相关性分析。
# cor.test
aa <- cor.test(e$a,e$b);aa
bb <- cor.test(d1$a1,d1$b1);bb输出结果：
当数据矩阵中存在大量缺失值时，可以利用corr.test()函数轻松化解难题。该函数可以直接对矩阵中的每两两变量进行相关性分析，同时能输出对应cor系数和p值。另外，还可以对等长度的单矩阵进行相关性分析，因此相比上面cor()和cor.test()函数，更加推荐corr.test()函数。但是大家一定要记住，该函数基于 psych R包，每次记得先library(psych)。废话不多说，进入例子学习：
Example 1 -- 单矩阵非等长度的相关性分析 # 单矩阵非等长度的相关性分析
?corr.test()
# corr.test
library(psych)
set.seed(3)
df1 <- rnorm(50, 1, 2)
df2 <- rnorm(100, 2, 3)
df3 <- runif(200, 10, 20)
df_all <- data.frame(df1, df2, df3);tail(df_all,20)
df <- df_all
# 需要将对应向量中非本身存在的数值替换成空值
df[51:200, 1] <- NA
df[101:200, 2] <- NA
tail(df,100)
str(df)
corr.test(df, method = "pearson")先构建了一个没有缺失的数据集，这里需要强调的是，当多个向量之间的长度不同时，如果将其用data.frame构建成数据框，系统会默认将短向量按当前向量数值的顺序补齐与其他向量长度保持一致。因此，如果我们要构建一个具有不等长的数据矩阵，需要对某列向量进行数值的缺失值替换，从而形成一个具有缺失值的数据框。 
从结果来看，1）没有报错；2）也符合我们的预期，数据矩阵中每两列之间都进行了相关性的比较；3）能看到样本量的大小，对应的cor系数和p值大小。因此，该函数得到的结果可读性比较高。 Example 2 -- 单矩阵等长度的相关性分析 （corr.test 函数）
这里就不写太多代码了，df_all这个数据集在上个例子中就已经展示了尾部5行，可以看到数据都是齐的，证明是等长度的数据矩阵。 # 等长度的数据矩阵相关性分析
str(df_all)
corr.test(df_all, method = "pearson")输出结果：首先看到样本数量大小为200，这是等长的。其次，能够详细的看到对应输出的结果，提高了整体的可读性。 
进一步，对得到的结果进行绘图
这里介绍两种专门绘制相关性矩阵分析的R包，1）GGally；2）corrplot。
# 绘图
#1 GGally 包
?ggcorr
library(GGally)
ggcorr(df,label=TRUE)
#2 corrplot包
?corrplot
library(corrplot)
res1 <- cor.mtest(df, conf.level=0.95,use="complete.obs")
N <- cor(df,use="complete.obs" ,method = "pearson") # use 是可以省略缺失值的相关性分析
corrplot(N,
p.mat=res1$p, # res1$p 中$表示提取符号 意为提取该数据的p列数据，p为显著性
insig ="label_sig", # 如果“label_sig”，标记与pch显著相关 用星号表示
sig.level=c(0.001, 0.01, 0.05), # 0.001 表示3个星号 以此类推
pch.col = "white",
pch.cex = 1) # 数字，颜色标签中的数字标签的cex关于结果具体完善代码就不细说了，有需要可以自行help然后查看描述及其例子。如果遇到问题，可以给推文评论或者私聊联系~
为什么说这是福利呢？因为从公司或者做实验得到的同类型但不同的数据矩阵很多，也就是有多个数据集。如何直接对两个数据矩阵之间做相关分析呢？目前好像没有函数可以直接做到，于是我的好兄弟和同学--谷大侠。他自己自定义一个能够计算两个矩阵等长度的相关性，原本计划是由他来写这篇推送的，但他最近忙于毕业，因此由我代劳。 # 自定义函数--计算两个矩阵之间的相关性
cor_function<-function(x,y,method){
library(reshape2)
if (is.null(y)) {
r <- cor(x, method = method)
sym <- TRUE
n <- t(!is.na(x)) %*% (!is.na(x))
t <- (r * sqrt(n - 2))/sqrt(1 - r^2)
p <- -2 * expm1(pt(abs(t), (n - 2), log.p = TRUE))
se <- sqrt((1 - r * r)/(n - 2))
nvar <- ncol(r)
p[p > 1] <- 1
r_re<-reshape2::melt(r)
r_re$value<-round(r_re$value,2)
p_re<-reshape2::melt(p)
Pz<-rep(1,nrow(p_re))
Pz[p_re$value<=0.001]<-'***'
Pz[p_re$value<=0.01&p_re$value>0.001]<-'**'
Pz[p_re$value<0.05&p_re$value>0.01]<-'*'
Pz[p_re$value>=0.05]<-NA
p_re<-cbind(p_re,Pz)
cor_data<-cbind(r_re,Pz)
}
else {
r <- cor(x, y, method = method)
sym = FALSE
n <- t(!is.na(x)) %*% (!is.na(y))
t <- (r * sqrt(n - 2))/sqrt(1 - r^2)
p <- -2 * expm1(pt(abs(t), (n - 2), log.p = TRUE))
se <- sqrt((1 - r * r)/(n - 2))
nvar <- ncol(r)
p[p > 1] <- 1
r_re<-reshape2::melt(r)
r_re$value<-round(r_re$value,2)
p_re<-reshape2::melt(p)
Pz<-rep(1,nrow(p_re))
Pz[p_re$value<=0.001]<-'***'
Pz[p_re$value<0.01&p_re$value>0.001]<-'**'
Pz[p_re$value<0.05&p_re$value>=0.01]<-'*'
Pz[p_re$value>=0.05]<-NA
p_re<-cbind(p_re,Pz)
cor_data<-cbind(r_re,Pz)
}
result <- list(cor_data=cor_data,r = r_re, n = n, t = t, p = p_re, se = se)
class(result) <- c("psych", "corr.test")
return(result)
}
# 构建两个矩阵
dat1 <- mtcars[,3:7]
dat2 <- iris[1:32,1:4]
str(dat1)
str(dat2)
# 画图 ----------------------------------------------------------------------
library(viridis)
library(ggplot2)
dat_cor <-cor_function(dat1,dat2,method='pearson')
ggplot(data=dat_cor$cor_data)+ # $cor_data 提取相关性r值
geom_tile(aes(x=Var1,y=Var2,fill=value))+ # x y 分别为dat_cor 中的Var1 Var2
scale_fill_viridis(option="viridis",limits=c(-1,1))+
theme(axis.text.x=element_text(hjust = 1,size=14),
axis.title.x=element_blank(),
axis.text.y=element_text(size=14),
axis.title.y=element_blank())+
geom_text(data=dat_cor$cor_data,aes(x=Var1,y=Var2,label=value,hjust=0.5))关于函数定义的部分就不展开讲解了，大家可以先copy进行使用。然后熟练掌握函数定义后就能理解了~ 结果输出：
|
