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