专注系列化、高质量的R语言教程
推文索引 | 联系小编 | 付费合集
我们最常用的相关系数是皮尔逊(Pearson)相关系数,也叫简单相关系数,用来衡量两个配对连续变量的线性相关程度。此外,还有斯皮尔曼(Spearman)相关系数和肯德尔(Kendall)相关系数可以度量有序变量之间的相关性。
在R语言中,这三个相关系数均可使用stats
工具包中的cor()
函数和cor.test()
函数进行计算和显著性检验。
cor(x,y=NULL,use="everything",method=c("pearson","kendall","spearman"))cor.test(x,y,alternative=c("two.sided","less","greater"),method=c("pearson","kendall","spearman"),exact=NULL,conf.level=0.95,continuity=FALSE,...)
本篇目录如下:
1 皮尔逊(Pearson)相关系数
1.1 计算公式
1.2 置换检验
1.3 贝塔检验
1.4 t检验
2 斯皮尔曼(Spearman)相关系数
2.1 秩
2.2 计算公式
1 皮尔逊(Pearson)相关系数
1.1 计算公式
使用、分别表示两个配对数值向量,则可以计算它们的协方差:
皮尔逊相关系数:
x<-c(44.4,45.9,41.9,53.3,44.7,44.1)y<-c(2.6,3.1,2.5,5.0,3.6,4.0)cor(x,y)##[1]0.826797
因为用来计算的数据只是来自总体的样本,其结果具有一定的“偶然性”,因此需要对其做显著性检验。
1.2 置换检验
相关系数的显著性检验的标准方式是使用置换试验。假设、的长度为,那么向量元素的可能配对方式共有(!表示阶乘)种(重复值不视为同一种情况),每种配对方式都可以计算出一个相关系数。使用表示这些相关系数中大于或等于实际样本计算的相关系数的个数,则显著性水平:
下面演示一下计算过程。
首先生成一个矩阵y2
,用来储存y
向量元素的所有可能排列情况。它的每一列表示一种排列情况,列数为。
library(e1071)all=permutations(6)y2=y[t(all)]y2=matrix(y2,nrow=6)
Note(s):这里使用的排列函数
permutations()
来自e1071
工具包。
然后再计算这720个相关系数并计算显著性水平:
all.r=cor(x,y2)sum(all.r>cor(x,y))/720##[1]0.04305556
Note(s):置换试验显示,相关系数的显著性水平等于0.043 < 0.05,因此可以有95%的把握认为两个向量具有正相关(即相关系数显著大于0)。
1.3 贝塔检验
使用置换试验进行显著性水平有明显的缺点,因为随着的增大,其计算量会迅速增大。有文章认为其只适合的情况[1]。
参考文献[2][2]指出,相关系数的概率密度函数可以使用下式表示:
其中,为贝塔函数,。
均值和方差:
则相关系数平方的累积概率为
进而的概率密度函数为
因此,服从贝塔分布(Beta Distribution):
在R中,可以使用pbeta()
函数计算贝塔分布的累积概率。在本例中,:
##显著性水平1-pbeta(cor(x,y)^2,0.5,2)##[1]0.04240094
1.4 t检验
t检验是对贝塔检验的近似。参考文献[2]指出,在z转换后近似服从自由度为的t分布。
因为,
所以,
z转换后,
在R语言中,可以使用pt()
函数计算t分布的累积概率:
t=cor(x,y)*sqrt(4)/sqrt(1-cor(x,y)^2)t##显著性水平(1-pt(t,4))*2##[1]0.04240094
cor.test()
函数就是使用t检验进行显著性检验的:
cor.test(x,y)##Pearson'sproduct-momentcorrelation####data:xandy##t=2.9397,df=4,p-value=0.0424##alternativehypothesis:truecorrelationisnotequalto0##95percentconfidenceinterval:##0.046308630.98046785##sampleestimates:##cor##0.826797
对于单变量线性回归来说,回归系数的显著性检验和相关系数的t检验的结果是一致的:
summary(lm(y~x))##Coefficients:##EstimateStd.ErrortvaluePr(>|t|)##(Intercept)-5.615113.09886-1.8120.1442##x0.198650.067582.9400.0424*
2 斯皮尔曼(Spearman)相关系数
对于有序分类变量来说,它的值不具有数值含义。虽然我们通常使用从1开始的正整数依次表示各个类别,但这些数值实际上只表示类别的顺序大小,而不能表示类别的差距大小,因此也可以使用大小顺序相同的其他任意数值表示。
2.1 秩
计算斯皮尔曼相关系数前需要先计算向量元素的“秩”(Rank)。对于一个向量,将它的元素按从小到达排序,则对于非重复元素来说,排序后的位置就是它的秩;对于重复元素来说,这些相等元素位置的平均数就是它们的秩。
例如,对下面的m
变量:
m=c(1,2,3,4,1,2,4)
将其元素按从小到大排列后:
sort(m)##[1]1122344
其中,两个1
的位置分别为1和2,则它们的秩均为1.5;两个4
的位置为6和7,则它们的秩为6.5;3的位置为5,且仅有一个3,则它的秩为5。
在R中,我们可以使用rank()
函数计算向量元素的秩:
rank(m,ties.method="average")##[1]1.53.55.06.51.53.56.5
2.2 计算公式
在得到向量元素的秩后,将秩视作普通数值,然后套用皮尔逊相关系数的公式就可以得到斯皮尔曼相关系数了:
m=c(1,2,3,4,1,2,4)n=c(1,3,4,2,1,3,3)mr=rank(m)nr=rank(n)##使用秩计算皮尔逊系数就是斯皮尔曼系数cor(mr,nr)##[1]0.5577955##使用原始数据计算斯皮尔曼系数cor(m,n,method="spearman")##[1]0.5577955
斯皮尔曼相关系数不仅仅只适用于有序分类变量。对于连续变量,若仅考虑其数值大小顺序,也可以计算它们的斯皮尔曼相关系数:
cor(x,y,method="spearman")##[1]0.6
当两个变量均不存在重复元素时,也可以使用下式计算斯皮尔曼相关系数:
由于上式中是秩差平方和,从而赋予了斯皮尔曼相关系数特定的含义,即它是使用配对元素在各自向量中排序的差别来度量相关性的。
m2=c(1,2,3,4,5,6,7)n2=c(1,3,4,6,7,2,5)cor(m2,n2,method="spearman")##[1]0.46428571-6*sum((m2-n2)^2)/(7^3-7)##[1]0.4642857
当向量中存在重复值时,上式计算结果就不同于斯皮尔曼系数了:
cor(m,n,method="spearman")##[1]0.55779551-6*sum((m-n)^2)/(7^3-7)##[1]0.8571429
由于斯皮尔曼相关系数的实质是使用秩计算的皮尔逊相关系数,因此前者的显著性检验也与后者类似:
cor.test(m,n,method="spearman")##Warningincor.test.default(m,n,method="spearman"):无法给连结计算精確p值####Spearman'srankcorrelationrho####data:mandn##S=24.763,p-value=0.1932##alternativehypothesis:truerhoisnotequalto0##sampleestimates:##rho##0.5577955cor.test(mr,nr,method="pearson")##Pearson'sproduct-momentcorrelation####data:mrandnr##t=1.5028,df=5,p-value=0.1932##alternativehypothesis:truecorrelationisnotequalto0##95percentconfidenceinterval:##-0.33669020.9231023##sampleestimates:##cor##0.5577955
参考资料
[1]
Algorithm AS 89: The Upper Tail Probabilities of Spearman's Rho:/10.2307/2347111
[2]
Significance Tests Which May be Applied to Samples from any Populations. II. The Correlation Coefficient Test:/10.2307/2983647