Pearson、Spearman、Kendall、Polychoric、Polyserial相關係數簡介及R計算

變量間Pearson、Spearman、Kendall、Polychoric、Tetrachoric、Polyserial、Biserial相關係數簡介及R計算對於給定數據集中，變量之間的關聯程度以及關係的方向，常通過相關係數衡量。就關係的強度而言，相關係數的值在+1和-1之間變化，值±1表示變量之間存在完美關聯程度，即完全相關時絕對值為1；隨著相關係數值趨於0，意味著變量之間的關係將減弱，完全不相關時為0。關係的方向由係數的符號表示；+號表示正向關係，-號表示負向關係。

圖示兩個變量之間的相關係數，正相關意味著圖表從左到右具有向上的斜率：隨著x值的增加，y值會變大；負相關性意味著圖表從左到右具有向下的斜率：隨著x值的增加，y值會變小；零（不相關）表示y不隨x的變化而變化。     

首先簡介常見的用於描述變量間相關性的係數，包括Pearson、Spearman、Kendall、Polychoric、Tetrachoric、Polyserial、(Point-)Biserial等。  

Pearson相關係數（皮爾森相關）是使用最廣泛的相關性統計量，用於測量兩組連續變量之間的線性關聯程度。

Pearson相關係數計算如下：

rxy，變量x和y的Pearson相關係數；

n，觀測對象的數量；

xi，x的第i個觀測值；

yi，y的第i個觀測值。

 

Pearson相關係數應用於連續變量，假定兩組變量均為正態分布、存在線性關係且等方差。線性關係假設兩個變量之間是線性響應的，等方差假設數據在回歸線上均勻分布。

  

Spearman秩相關係數（斯皮爾曼等級相關）是一種非參數統計量，其值與兩組相關變量的具體值無關，而僅僅與其值之間的大小關係有關。Spearman秩相關依據兩列成對等級的各對等級數之差進行計算，所以又稱為「等級差數法」。當變量在至少是有序的尺度上測量時，它是合適的相關分析方法。

Spearman秩相關係數計算如下：

ρ，Spearman秩相關係數；

di，對應變量的秩之差，即兩個變量分別排序後成對的變量位置（等級）差；

n，觀測對象的數量。

 

Spearman秩相關同樣應用於連續變量，與Pearson相關相比Spearman秩相關不要求變量的正態性和等方差假設，且對異常值的敏感度較低（該方法基於變量的排序，因此異常值的秩次通常不會有明顯變化），因此適用範圍通常更廣。但方法較為保守，統計效能較Pearson相關係數低，容易忽略一些不太強的線性關係。

此外，Spearman秩相關要求數據必須至少是有序的，一個變量的得分必須與另一個變量單調相關（monotonically related）。 

 

Kendall 相關係數則用於計算分類變量間的秩相關，用於反映分類變量相關性的指標，適用於兩個分類變量均為有序分類的情況。

考慮兩組變量，x和y，它們各自的觀測值數量均為n，則x與y觀測值可能配對的總數為n(n-1)/ 2。由於x和y為分類變量，需要首先根據類別表示的重要度人工賦值。隨後考察x和y的關係對，如果xi<yi且xj<yj，或xi>yi且xj>yj，則該關係對是一致的（concordant），反正則不一致（discordant）。一致關係對數量與不一致關係對數量的差值除以總關係對數量，可得Kendall 相關係數：

如果一致對的數量比不一致對的數量大得多，則變量是正相關的；如果一致對的數目比不一致對的數目少得多；則變量是負相關的；如果一致對的數目與不一致對的數目大致相同，則變量之間的關係很弱。

  

Polychoric相關（多分格相關）度量多個對象之間關於有序變量（有時稱為「有序類別」數據）之間的一致性。當以列聯表的形式組織數據時，兩個分類自變量被排序，據此計算Polychoric相關係數。

對於2×2列聯表的情況，Polychoric相關係數也稱為Tetrachoric相關係數（作為Polychoric相關的一種常見類型）。通過以下對Tetrachoric相關的描述即可理解Polychoric相關的定義。  

Tetrachoric相關（二元有序變量間的相關，Polychoric相關的某種常見類型）

Tetrachoric相關（四分相關）是在二元正態性假設下從2×2表推斷出的Pearson相關，用於測量二元數據一致性。Tetrachoric相關要求基本變量來自正態分布，並且二元數據中存在一個潛在的連續梯度，即觀測值的特徵應該是連續而非離散的。

首先將觀察數值矩陣獲得列聯表，並通過下式計算：

 

 

Polyserial相關（定量變量和序數變量的相關）

Polyserial相關（多系列相關）測量的是兩個連續變量之間的相關關係，它們具有二元正態分布，其中一個變量可以直接觀測到（以定量數值記錄），而另一個變量無法被觀測（以序數值記錄）。通過將可觀測的連續變量分類為有限的離散有序值集，可以從可觀測的有序變量獲得不可觀測有序變量的信息。

Biserial相關（連續變量和二元有序變量的相關，Polyserial相關的某種特例）

Biserial相關係數為Polyserial相關的一種特例，用於測量一組連續變量和一組二元變量的線性關係，二元變量是二分序數類型，具有潛在的連續性。

Y0，x=0時變量對的平均分；

Y1，x=1時變量對的平均分；

p，x=1時變量對的比例；

q，x=0時變量對的比例；

σy，總體標準偏差。

  

Point-Biserial相關（連續變量和二元分類變量的相關）

與Biserial相關係數相比，Point-Biserial相關係數用於測量一組連續變量和一組二元分類變量的線性關係，分類變量是無序的。

M1，二元變量組「1」對象對應的連續變量的均值；

M0，二元變量組「0」對象對應的連續變量的均值；

Sn，連續變量的標準偏差；

p，二元變量組「1」對象所佔總對象的比例；

q，二元變量組「0」對象所佔總對象的比例。

接下來展示在R中計算上述提到的相關係數的方法。  

Pearson、Spearman和Kendall相關

在R中，cor()可用於計算Pearson、Spearman和Kendall相關矩陣，cov()可用於獲得協方差矩陣。

##Pearson、Spearman、Kendall 相關
data(mtcars)

#標準化不影響相關係數計算值，但可以讓數據服從均值 0，標準差 1 的等方差結構
mtcars <- scale(mtcars)

#協方差計算，cov()
cov_pearson <- cov(mtcars, method = 'pearson')
cov_pearson

cov_spearman <- cov(mtcars, method = 'spearman')
cov_spearman

cov_kendall <- cov(mtcars, method = 'kendall')
cov_kendall

#相關係數計算，cor()
cor_pearson <- cor(mtcars, method = 'pearson')
cor_pearson

cor_spearman <- cor(mtcars, method = 'spearman')
cor_spearman

cor_kendall <- cor(mtcars, method = 'kendall')
cor_kendall

#相關圖，例如
library(corrplot)

corrplot(cor_pearson, method = 'number', number.cex = 0.8, diag = FALSE, tl.cex = 0.8)
corrplot(cor_pearson, add = TRUE, type = 'upper', method = 'pie', diag = FALSE, tl.pos = 'n', cl.pos = 'n')

#輸出，例如
write.table(cor_pearson, 'cor_pearson.txt', sep = '\t', col.names = NA, quote = FALSE)

直接指定數據集，默認計算所有變量間的相關係數，獲得斜對角線對稱的矩陣。

也可指定兩組變量集，獲得相互之間兩兩變量間非對稱的相關矩陣。

#指定兩組變量集，獲得非對稱的相關矩陣，例如
x <- c('mpg', 'cyl', 'disp', 'hp')
y <- c('drat', 'wt', 'qsec')

cor_pearson_xy <- cor(mtcars[x], mtcars[y], method = 'pearson')
cor_pearson_xy

#相關圖
corrplot(cor_pearson_xy, method = 'square', addCoef.col = 'black', number.cex = 0.8, tl.cex = 1.2)

#輸出，例如
write.table(cor_pearson_xy, 'cor_pearson_xy.txt', sep = '\t', col.names = NA, quote = FALSE)

  

偏相關是指在控制一個或多個定量變量時，另外兩個定量變量之間的相互關係。R包ggm中提供的命令pcor()可以計算偏相關係數。

##偏相關，ggm 包 pcor()
library(ggm)

#要計算相關係數的兩個變量，或指定下標
x1 <- c('mpg', 'cyl')

#要控制的條件變量，或指定下標
x2 <- c('drat', 'wt', 'qsec')

#指定協方差矩陣，計算偏相關
pcor_pearson <- pcor(c(x1, x2), cov(mtcars, method = 'pearson'))
pcor_pearson

  

psych包提供了計算這些相關係數的方法。

psych包也能計算Polyserial和Biserial相關，但文檔中沒提供示例，沒看明白……

##Polychoric、Tetrachoric
library(psych)

#Polychoric 相關
data(bock)

polyc <- polychoric(lsat6)
polyc

#Tetrachoric 相關
tetr <- tetrachoric(lsat6[ ,1:2])
tetr

  

Polyserial和(Point-)Biserial相關

以ltm包提供的方法為例。

##Polyserial、(Point-)Biserial
library(ltm)

#Polyserial 相關
mpg <- subset(ggplot2::mpg, class == 'midsize' | class == 'compact')
polys <- polyserial(mpg$cty, mpg$class, std.err = TRUE)
polys

#Point-Biserial 相關
poi_biser <- biserial.cor(mpg$cty, mpg$class)
poi_biser

#Biserial 相關
biser <- biserial.cor(mtcars$mpg, mtcars$vs)
biser

通常來講，相關性分析是一種用於描述變量關聯程度的探索性分析方法，而非確立因果關係的模型，不涉及假設檢驗過程。但如果有必要，仍可以計算相關係數的顯著性，評估哪些變量間的關聯程度是更重要的。

一些R包提供了計算變量間相關係數顯著性的方法。此外，也可以自寫函數獲得，見下文。  

計算相關矩陣及顯著性水平。

library(psych)

#所有變量間相關係數的對稱矩陣
corr_matrix <- corr.test(mtcars, method = 'pearson')
corr_matrix$r #相關矩陣
corr_matrix$p #p 值矩陣

#相關圖，只展示 p < 0.05 的相關係數
library(corrplot)

col1 <- colorRampPalette(c('blue4', 'blue', 'white', 'orange', 'red3'))
corrplot(corr_matrix$r, p.mat = corr_matrix$p, sig.level = 0.05, insig = 'blank', method = 'number',
diag = FALSE, col = col1(21), tl.cex = 1)
corrplot(corr_matrix$r, p.mat = corr_matrix$p, sig.level = 0.05, insig = 'blank', method = 'circle',
add = TRUE, type = 'upper', diag = FALSE, col = col1(21), tl.pos = 'n', cl.pos = 'n')

#自定義篩選，例如選擇 |r| >=0.7，p < 0.05 的結果，將不滿足條件的相關係數值賦值為 0 後輸出
corr_matrix$p[corr_matrix$p >= 0.05] <- -1
corr_matrix$p[corr_matrix$p < 0.05 & corr_matrix$p >= 0] <- 1
corr_matrix$p[corr_matrix$p == -1] <- 0

corr_matrix$r[abs(corr_matrix$r) < 0.7] <- 0
corr_matrix$r <- corr_matrix$r * corr_matrix$p
write.table(corr_matrix$r, 'corr_matrix_select.txt', sep = '\t', col.names = NA, quote = FALSE)

#給定兩組變量間相關係數的非對稱矩陣
x <- as.matrix(mtcars[c('mpg', 'cyl', 'disp', 'hp')])
y <- as.matrix((mtcars[c('drat', 'wt', 'qsec')])

corr_matrix <- corr.test(x, y, method = 'pearson')
corr_matrix$r #相關矩陣
corr_matrix$p #p 值矩陣

#相關圖，只展示 p < 0.05 的相關係數
col1 <- colorRampPalette(c('blue4', 'blue', 'white', 'orange', 'red3'))
corrplot(corr_matrix$r, p.mat = corr_matrix$p, sig.level = 0.05, insig = 'blank',
method = 'square', addCoef.col = 'black', col = col1(21), number.cex = 0.8, tl.cex = 1.2)

  

計算相關矩陣及顯著性水平。

library(Hmisc)

#所有變量間相關係數的對稱矩陣
rcorr_matrix <- rcorr(as.matrix(mtcars), type = 'pearson')
rcorr_matrix$r #相關矩陣
rcorr_matrix$P #p 值矩陣

#給定兩組變量間相關係數的非對稱矩陣
x <- as.matrix(mtcars[c('mpg', 'cyl', 'disp', 'hp')])
y <- as.matrix((mtcars[c('drat', 'wt', 'qsec')]))

rcorr_matrix <- rcorr(x, y, type = 'pearson')
rcorr_matrix$r #相關矩陣
rcorr_matrix$P #p 值矩陣

#相關圖、自定義結果篩選等，參考上述

  

置換檢驗是個百搭的非參數檢驗方法，相關係數的顯著性可根據置換檢驗的原理獲得。

上述提到的所有相關係數，包括Polychoric、Tetrachoric、Polyserial、(Point-)Biserial等，如果找不到計算顯著性的R包，不妨考慮手寫函數計算，其實並不難。

#計算觀測值的相關係數（cor0），還是以 Pearson 相關為例，其它類似
cor0 <- cor(mtcars, method = 'pearson')

#隨機置換數據 999 次，計算每次置換後數據計算的相關係數（corN），並統計 |corN|>|cor0| 的頻數
p_num <- cor0
p_num[abs(p_num)>0] <- 1

set.seed(123)
for (i in 1:999) {
random <- apply(mtcars, 2, sample)
corN <- cor(random, method = 'pearson')

corN[abs(corN) >= abs(cor0)] <- 1
corN[abs(corN) < abs(cor0)] <- 0
p_num <- p_num + corN
}

#p 值矩陣，即 |corN|>|cor0| 的概率
p <- p_num/1000
p

#相關圖比較，僅顯著（p < 0.05）的相關係數標以背景色
#左圖為手寫的置換檢驗結果，右圖為 psych 包獲得的結果，二者是一致的
library(corrplot)
library(psych)
cor_psych <- corr.test(mtcars, method = 'pearson')

layout(matrix(c(1,2), 1, 2, byrow = TRUE))
corrplot(cor0, method = 'square', type = 'lower', p.mat = p, sig.level = 0.05, insig = 'blank',
addCoef.col = 'black', diag = FALSE, number.cex = 0.8, tl.cex = 0.8)
corrplot(cor_psych$r, method = 'square', type = 'lower', p.mat = cor_psych$p, sig.level = 0.05, insig = 'blank',
addCoef.col = 'black', diag = FALSE, number.cex = 0.8, tl.cex = 0.8)

Pearson, Spearman & Kendall：https://www.statisticssolutions.com/correlation-pearson-kendall-spearman/

(Point-)Biserial：https://www.statisticshowto.datasciencecentral.com/point-biserial-correlation/

Tetrachoric & Polychoric：http://john-uebersax.com/stat/tetra.htm

Polyserial：http://support.sas.com/documentation/cdl/en/procstat/63963/HTML/default/viewer.htm#procstat_corr_sect019.htm