https://corrr.tidymodels.org/articles/using-corrr.html
https://easystats.github.io/correlation/
相关性:判断两个变量之间的相关关系强度和方向,无论是否独立。
分类变量
如果独立性检验的结果表明两个变量之间不独立,那么如何量化它们之间相关性的强弱?
Phi 系数、列联系数和 Cramer’s V 系数
vcd 包里的函数 assocstats( )
可以用来计算列联表的 Phi 系数、列联系数和 Cramer’s V 系数。其中, Phi 系数只适用于四格表。
Code
library(vcd)
mytable <- table(Arthritis$Sex, Arthritis$Treatment)
assocstats(mytable)
#> X^2 df P(> X^2)
#> Likelihood Ratio 0.73748 1 0.39047
#> Pearson 0.73653 1 0.39078
#>
#> Phi-Coefficient : 0.094
#> Contingency Coeff.: 0.093
#> Cramer's V : 0.094
Kappa 统计量
对于配对列联表,可以计算一致性指标 Kappa 统计量。 epiDisplay 包里的函数 kap( )可以用于计算一致性的比例以及 Kappa 统计量的值
Code
my.matrix <- matrix(c(11, 2, 12, 33), nrow = 2)
sum(my.matrix)
#> [1] 58
vcd::Kappa(my.matrix)
#> value ASE z Pr(>|z|)
#> Unweighted 0.455 0.1153 3.945 7.97e-05
#> Weighted 0.455 0.1153 3.945 7.97e-05
epiDisplay::kap(my.matrix)
#>
#> Table for calculation of kappa
#> A B
#> A 11 12
#> B 2 33
#>
#> Observed agreement = 75.86 %
#> Expected agreement = 55.71 %
#> Kappa = 0.455
#> Standard error = 0.121 , Z = 3.762 , P value = < 0.001
#>
共 58 个对象,每一对象用两种检测方法检测,其中11 个对象的两种检测结果都为阳性, 33 个对象的两种检测结果都是阴性,所以总一致性为 (11 + 33)/58 ≈ 75.86% 。
Code
chisq.test(my.matrix)$expected
#> [,1] [,2]
#> [1,] 5.155172 17.84483
#> [2,] 7.844828 27.15517
为了解释期望一致性和 Kappa 值的含义,先计算各个单元格的期望频数。 对角线上的这两个单元格对应的期望频数分别约为5.155172 和27.15517 ,因此期望一致性为 (5.155172+27.15517)/58≈ 55.71% 。期望一致性是假定两种方法的检测结果都是完全随机的情况下的 一致性。也就是说,即使两种检测方法都毫无作用,平均也能达到 55.71% 的一致性。 Kappa 统计量是超出随机的一致性的部分占最大可能超出随机的一致性的比例。在本例中,前者为 75.86% − 55.71% , 后者为 100% − 55.71% 。 因此, Kappa 值为 (75.86 - 55.71)/(100 - 55.71) ≈ 0.455
马赛克图
马赛克图中的矩形面积正比于多维列联表中单元格的频率
Code
mosaicplot(mytable,xlab ="Sex",ylab = "Treatment",las = 1)
Code
连续变量
如果两个连续变量不相互独立时,使用协方差(covariance)来描述两个变量的关系。
协方差(或相关系数)为零,不相关,不存在线性关系,但可能存在非线性关系。
Code
df <- mpg[,c(3,8,9)]
cov(df) # 协方差矩阵
#> displ cty hwy
#> displ 1.669158 -4.39069 -5.893111
#> cty -4.390690 18.11307 24.225432
#> hwy -5.893111 24.22543 35.457779
相关系数
相关系数的取值范围: \([-1,1]\)
\[
r(X,Y)=\frac {\sum_{i=1}^n (x_i-\bar x)(y_i-\bar y)}{\sqrt{\sum_{i=1}^n (x_i-\bar x)^2 \sum_{i=1}^n (y_i-\bar y)^2}}
\]
Code
# Pearson's 积差相关系数 一般要求两个连续变量都服从正态分布
cor(df,use = "everything",method="pearson") # default
#> displ cty hwy
#> displ 1.000000 -0.7985240 -0.7660200
#> cty -0.798524 1.0000000 0.9559159
#> hwy -0.766020 0.9559159 1.0000000
correlation::correlation(df,method = "pearson",p_adjust = "holm")
displ |
cty |
-0.7985240 |
0.95 |
-0.8406782 |
-0.7467508 |
-20.20515 |
232 |
0 |
Pearson correlation |
234 |
displ |
hwy |
-0.7660200 |
0.95 |
-0.8142727 |
-0.7072539 |
-18.15085 |
232 |
0 |
Pearson correlation |
234 |
cty |
hwy |
0.9559159 |
0.95 |
0.9433129 |
0.9657663 |
49.58470 |
232 |
0 |
Pearson correlation |
234 |
Code
# Spearman's rank相关系数 非参数
cor(df,method = "spearman")
#> displ cty hwy
#> displ 1.0000000 -0.8809049 -0.8266576
#> cty -0.8809049 1.0000000 0.9542104
#> hwy -0.8266576 0.9542104 1.0000000
correlation::correlation(df,method = "spearman",p_adjust = "holm")
displ |
cty |
-0.8809049 |
0.95 |
-0.9073926 |
-0.8474471 |
4016568.96 |
0 |
Spearman correlation |
234 |
displ |
hwy |
-0.8266576 |
0.95 |
-0.8643401 |
-0.7797446 |
3900726.77 |
0 |
Spearman correlation |
234 |
cty |
hwy |
0.9542104 |
0.95 |
0.9406973 |
0.9647003 |
97781.21 |
0 |
Spearman correlation |
234 |
Code
# Kendall's tau相关系数 非参数
cor(df,method = "kendall")
#> displ cty hwy
#> displ 1.0000000 -0.7210828 -0.6536974
#> cty -0.7210828 1.0000000 0.8628045
#> hwy -0.6536974 0.8628045 1.0000000
correlation::correlation(df,method = "kendall",p_adjust = "holm")
displ |
cty |
-0.7210828 |
0.95 |
-0.7596258 |
-0.6774923 |
-15.53533 |
0 |
Kendall correlation |
234 |
displ |
hwy |
-0.6536974 |
0.95 |
-0.6999287 |
-0.6020109 |
-14.13933 |
0 |
Kendall correlation |
234 |
cty |
hwy |
0.8628045 |
0.95 |
0.8392948 |
0.8830936 |
18.39871 |
0 |
Kendall correlation |
234 |
相关图(correlogram)
Code
ggcorrplot::ggcorrplot(
corr = cor(df,use = "everything",method="pearson") ,
lab = T
)
显著性检验
零假设为变量之间不相关(即两个总体的相关系数为 0 ) 。函数 cor.test( )
可用于对相关系数进行显著性检 验。
统计量
\[
t=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}
\]
Code
cor.test(df$displ,df$hwy)
#>
#> Pearson's product-moment correlation
#>
#> data: df$displ and df$hwy
#> t = -18.151, df = 232, p-value < 2.2e-16
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#> -0.8142727 -0.7072539
#> sample estimates:
#> cor
#> -0.76602
psych包corr.test()
计算相关系数矩阵和显著性检验
Code
psych::corr.test(df)
#> Call:psych::corr.test(x = df)
#> Correlation matrix
#> displ cty hwy
#> displ 1.00 -0.80 -0.77
#> cty -0.80 1.00 0.96
#> hwy -0.77 0.96 1.00
#> Sample Size
#> [1] 234
#> Probability values (Entries above the diagonal are adjusted for multiple tests.)
#> displ cty hwy
#> displ 0 0 0
#> cty 0 0 0
#> hwy 0 0 0
#>
#> To see confidence intervals of the correlations, print with the short=FALSE option
print(psych::corr.test(df), short = FALSE)
#> Call:psych::corr.test(x = df)
#> Correlation matrix
#> displ cty hwy
#> displ 1.00 -0.80 -0.77
#> cty -0.80 1.00 0.96
#> hwy -0.77 0.96 1.00
#> Sample Size
#> [1] 234
#> Probability values (Entries above the diagonal are adjusted for multiple tests.)
#> displ cty hwy
#> displ 0 0 0
#> cty 0 0 0
#> hwy 0 0 0
#>
#> Confidence intervals based upon normal theory. To get bootstrapped values, try cor.ci
#> raw.lower raw.r raw.upper raw.p lower.adj upper.adj
#> displ-cty -0.84 -0.80 -0.75 0 -0.85 -0.74
#> displ-hwy -0.81 -0.77 -0.71 0 -0.81 -0.71
#> cty-hwy 0.94 0.96 0.97 0 0.94 0.97