|
相關(guān)性分析通常用來定量描述兩個變量之間的聯(lián)系，正相關(guān)？負(fù)相關(guān)？不存在相關(guān)？等。常見相關(guān)性指標(biāo)：Pearson相關(guān)系數(shù)（參數(shù)檢驗），Spearman及Kendall（非參數(shù)檢驗）。本期與分享R語言相關(guān)分析時的小技巧：1) 一個矩陣中等長度兩兩變量的相關(guān)性，常用到cor() 函數(shù)和cor.test() 函數(shù)；2) 一個矩陣中非等長度兩個變量的相關(guān)性，常用corr.test()函數(shù)；3) 兩個矩陣中變量對應(yīng)的相關(guān)性，這里編了一個函數(shù)與大家分享。 (本期不做偏相關(guān)分析的講解）1）Pearson相關(guān)系數(shù)（積差相關(guān)系數(shù)）通常用來表示相關(guān)性大小中最常用的指標(biāo)。cor系數(shù)介于-1~1之間，越接近0相關(guān)性越弱，越接近-1 or 1則相關(guān)性越強，并存在對應(yīng)的負(fù)正相關(guān)。2）Spearman等級相關(guān)系數(shù)（秩相關(guān)系數(shù)）利用兩個變量的秩次大小進(jìn)行比較，屬于非參數(shù)檢驗。3）Kendall等級相關(guān)系數(shù)（和諧系數(shù)）也屬于一種非參數(shù)檢驗。等長度是兩個變量的樣本量相等。主要利用 cor() 函數(shù) 和 cor.test() 函數(shù)對不同方法的講解。以一組滿足正態(tài)分布的連續(xù)變量和不滿足正態(tài)分布的有序變量進(jìn)行對比。例子講解前，我們先用？cor（）/ help（cor）了解一下函數(shù)的用法和命令。
函數(shù)的基本命令：cor(x, use= , method= )x為一個數(shù)值向量，矩陣或者數(shù)據(jù)框都可以； y這里默認(rèn)為NULL； na.rm 表示缺失值的邏輯，計算比較時是否保留缺失值（T表示去除缺失值，但是這里只能用于var函數(shù)中）；
use 指定缺失數(shù)據(jù)的處理方式。可選項：all.obs(假設(shè)不存在缺失數(shù)據(jù))、everything(數(shù)據(jù)存在缺失值時，相關(guān)系數(shù)計算結(jié)果會顯示missing)、complete.obs(行刪除)、pairwise.complete.obs(成對刪除)； method 指定相關(guān)系數(shù)的類型?？蛇x類型為pearson、spearman、kendall。
清楚了函數(shù)對應(yīng)命令的參數(shù)后，進(jìn)行例子的理解 Example 1 -- 滿足正態(tài)分布的連續(xù)變量利用str（）檢查數(shù)據(jù)結(jié)構(gòu)的結(jié)果，num表示數(shù)值型的連續(xù)變量。#1 構(gòu)建一組符合正態(tài)分布的數(shù)據(jù)
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)
# 將上述構(gòu)建的數(shù)據(jù)組合成數(shù)據(jù)框, 并檢查數(shù)據(jù)結(jié)構(gòu)（num數(shù)值型）
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輸出結(jié)果：當(dāng)使用參數(shù)與非參數(shù)的方法對符合正態(tài)分布的連續(xù)變量進(jìn)行相關(guān)性比較時，結(jié)果是存在差異的。Pearson為一個結(jié)果，而Spearman和kendall為兩一個結(jié)果。
Example 2 -- 不滿足正態(tài)分布的有序變量 利用as.integer（）將向量轉(zhuǎn)換成整數(shù)型的，使其變成有序變量，數(shù)據(jù)結(jié)構(gòu)為int，表示整數(shù)型。 #2 構(gòu)建一組不滿足正態(tài)分布的數(shù)據(jù)
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輸出結(jié)果：從結(jié)果也能發(fā)現(xiàn)不同方法進(jìn)行的相關(guān)性分析確實存在差異，但是具體哪種更好，可能需要大家多做嘗試。 
接下來，利用cor.test（）函數(shù)對上面例子做詳細(xì)展示，可以看得到對應(yīng)的p值大小，有利于大家解讀和后續(xù)可視化。還是先咨詢下該函數(shù)的用法和對應(yīng)參數(shù)。 
函數(shù)的基本命令：cor.test(x, y, alternative= , method= ,...) 
這里不做詳細(xì)介紹，有問題的話可以評論留言~但這里需要注意的是：不能直接用數(shù)據(jù)框進(jìn)行分析，而是需要用兩兩向量，因此用$提取對應(yīng)向量進(jìn)行相關(guān)性分析。
# cor.test
aa <- cor.test(e$a,e$b);aa
bb <- cor.test(d1$a1,d1$b1);bb輸出結(jié)果：
當(dāng)數(shù)據(jù)矩陣中存在大量缺失值時，可以利用corr.test()函數(shù)輕松化解難題。該函數(shù)可以直接對矩陣中的每兩兩變量進(jìn)行相關(guān)性分析，同時能輸出對應(yīng)cor系數(shù)和p值。另外，還可以對等長度的單矩陣進(jìn)行相關(guān)性分析，因此相比上面cor()和cor.test()函數(shù)，更加推薦corr.test()函數(shù)。但是大家一定要記住，該函數(shù)基于 psych R包，每次記得先library(psych)。廢話不多說，進(jìn)入例子學(xué)習(xí)：
Example 1 -- 單矩陣非等長度的相關(guān)性分析 # 單矩陣非等長度的相關(guān)性分析
?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
# 需要將對應(yīng)向量中非本身存在的數(shù)值替換成空值
df[51:200, 1] <- NA
df[101:200, 2] <- NA
tail(df,100)
str(df)
corr.test(df, method = "pearson")先構(gòu)建了一個沒有缺失的數(shù)據(jù)集，這里需要強調(diào)的是，當(dāng)多個向量之間的長度不同時，如果將其用data.frame構(gòu)建成數(shù)據(jù)框，系統(tǒng)會默認(rèn)將短向量按當(dāng)前向量數(shù)值的順序補齊與其他向量長度保持一致。因此，如果我們要構(gòu)建一個具有不等長的數(shù)據(jù)矩陣，需要對某列向量進(jìn)行數(shù)值的缺失值替換，從而形成一個具有缺失值的數(shù)據(jù)框。 
從結(jié)果來看，1）沒有報錯；2）也符合我們的預(yù)期，數(shù)據(jù)矩陣中每兩列之間都進(jìn)行了相關(guān)性的比較；3）能看到樣本量的大小，對應(yīng)的cor系數(shù)和p值大小。因此，該函數(shù)得到的結(jié)果可讀性比較高。 Example 2 -- 單矩陣等長度的相關(guān)性分析 （corr.test 函數(shù)）
這里就不寫太多代碼了，df_all這個數(shù)據(jù)集在上個例子中就已經(jīng)展示了尾部5行，可以看到數(shù)據(jù)都是齊的，證明是等長度的數(shù)據(jù)矩陣。 # 等長度的數(shù)據(jù)矩陣相關(guān)性分析
str(df_all)
corr.test(df_all, method = "pearson")輸出結(jié)果：首先看到樣本數(shù)量大小為200，這是等長的。其次，能夠詳細(xì)的看到對應(yīng)輸出的結(jié)果，提高了整體的可讀性。 
進(jìn)一步，對得到的結(jié)果進(jìn)行繪圖
這里介紹兩種專門繪制相關(guān)性矩陣分析的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 是可以省略缺失值的相關(guān)性分析
corrplot(N,
p.mat=res1$p, # res1$p 中$表示提取符號 意為提取該數(shù)據(jù)的p列數(shù)據(jù)，p為顯著性
insig ="label_sig", # 如果“l(fā)abel_sig”，標(biāo)記與pch顯著相關(guān) 用星號表示
sig.level=c(0.001, 0.01, 0.05), # 0.001 表示3個星號 以此類推
pch.col = "white",
pch.cex = 1) # 數(shù)字，顏色標(biāo)簽中的數(shù)字標(biāo)簽的cex關(guān)于結(jié)果具體完善代碼就不細(xì)說了，有需要可以自行help然后查看描述及其例子。如果遇到問題，可以給推文評論或者私聊聯(lián)系~
福利大贈送~自定義函數(shù)計算兩個矩陣等長度的相關(guān)性 為什么說這是福利呢？因為從公司或者做實驗得到的同類型但不同的數(shù)據(jù)矩陣很多，也就是有多個數(shù)據(jù)集。如何直接對兩個數(shù)據(jù)矩陣之間做相關(guān)分析呢？目前好像沒有函數(shù)可以直接做到，于是我的好兄弟和同學(xué)--谷大俠。他自己自定義一個能夠計算兩個矩陣等長度的相關(guān)性，原本計劃是由他來寫這篇推送的，但他最近忙于畢業(yè)，因此由我代勞。 # 自定義函數(shù)--計算兩個矩陣之間的相關(guān)性
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)
}
# 構(gòu)建兩個矩陣
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 提取相關(guān)性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))關(guān)于函數(shù)定義的部分就不展開講解了，大家可以先copy進(jìn)行使用。然后熟練掌握函數(shù)定義后就能理解了~ 結(jié)果輸出：
|
