1 Pearson 相关系数\(r\)定义
在统计学里，两个序列（变量）的相关系数是Pearson 相关系数\(r\)，或Spearman相关系数\(\rho\)，或Kendall相关系数\(\tau\) 。 期中两个连续变量的Pearson 相关系数\(r\)使用最多，其定义是：
\[ r=\frac{Cov(X,Y)}{\sqrt{Var(X)\cdot Var(Y)}} \tag{1}\]
公式 Equation 1 说明相关系数是将协方差对标准差（方差的平方根）进行归一化的数值。 下图是变量X, Y样本\((x_i, y_i|i=1,2,...,20)\)的散点图，线性回归结果，Pearson r值及其显著性水平（p值）。
2 关于\(r\)的常见错误理解和解释
2.1 \(r=0\) 是否表明两个变量不相关？ No
Pearson 相关系数只描述变量的线性相关性。\(r=0\)或者r值很小，只能说明两个变量没有或者具有很弱的线性相关性，但不能排除两者是否具有非线性关系。 下图很直观地说明这个问题：
Code
set.seed(10)
x1=seq(from=1,to=20,by=1)
y1=2*x1+rnorm(20,0,10)
x2=seq(from=1,to=20,by=1)
y2=(x2-10.5)^2+rnorm(20,0,10)
par(mfrow=c(1,2))
plot(x1,y1,pch=1)
mod1 <- lm(y1 ~ x1)
cf <- coef(mod1)
eq1 = paste0("y1 = ", round(cf[2],1), "*x1 ", round(cf[1],1))
title(main = eq1,col.main="blue")
abline(mod1, col="blue")
pcc1=paste0("Pearson r = ",round(cor(x1,y1),2),"\n p=",format(cor.test(x1,y1)$p.value,scientific = TRUE,digits=3))
text(15, 2, pcc1, cex = 1,col="blue")
plot(x2,y2,pch=2)
mod2 <- lm(y2 ~ x2)
cf <- coef(mod2)
eq2 = paste0("y2 = ", round(cf[2],1), "*x2+", round(cf[1],1))
title(main = eq2,col.main="blue")
abline(mod2, col="blue")
pcc2=paste0("Pearson r = ",round(cor(x2,y2),2),"\n p=",format(cor.test(x2,y2)$p.value,scientific = TRUE,digits=3))
text(10, 60, pcc2, cex = 1,col="blue")
上图右图对应的样本来自于\(y_i=(x_i-10.5)^2+\epsilon_i\) , \(\epsilon \sim N(0,\sqrt{10})\)，随机变量X, Y 存在二次函数关系，但线性关系很弱。
结论：当\(r\) 为零，或者非常小的时候，我们只能说变量之间没有或者存在非常小的线性相关性。
2.2 r不等于零，比如r>0.5是否表示变量之间存在相关性？ No
经常有人在获得一个中等强度以上的Pearson r值，并且具有显著性（比如p<0.05）, 便主张变量之间存在线性关系。
前面说过，Pearson r描述变量间的线性相关性，并不表示变量之间是线性关系。 因为即使r值很大的情况下（比如r>0.8）,也有可能实际关系并不是线性关系，如下图：
Code
set.seed(10)
x=seq(from=1,to=20,by=1)
y=x^2+rnorm(20,0,10)
# plotting the regression line
mod <- lm(y ~ x)
cf <- coef(mod)
plot(x,y)
eq = paste0("y = ", round(cf[2],1), "*x ", round(cf[1],1))
title(main = eq,col.main="blue")
abline(mod, col="blue")
pcc=paste0("Pearson r = ",round(cor(x,y),2),"\n p=",format(cor.test(x,y)$p.value,scientific = TRUE,digits=3))
text(16, 50, pcc, cex = 1,col="blue")
text(5, 250, "The true relation is : y=x^2", cex = 1,col="black")
2.3 相关系数的显著性水平p<0.05， 或者p<0.01 (比如：r=0.35, p=0.012), 能说明变量之间存在很强的相关性吗？ No.
显著性水平p值和相关性强弱是两个不同的概念， 相关性强弱体现在r值上，通常的经验值是（thumb of rules）:
\(r<0.1\) 非常弱线性相关
\(0.1<r<0.4\) 弱线性相关性
\(0.4<r<0.6\) 中等线性相关性
\(0.6<r<0.8\) 强线性相关性
\(r>0.8\) 极强线性相关性
而p值是对r值的显著性水平的检验， 是一个单样本的t-检验，这个显著性检验的零假设和备择假设分别是：
\(H_0: r=0\)
\(H_1: r\ne 0\)
所以，你如果得到一个r值的显著性水平\(p<0.05\)，按照零假设检验的方法，你可以拒绝零假设 \(H_0\) , 接受备择假设\(H_1\)。 比如\(r=0.8, p=0.01\) , 并不能得到"存在显著的强相关性"的结论， 因为这个"显著性检验"并没有检验相关性的"强弱"， 拒绝\(H_0\), 你接受的\(H_1\) 是存在 \(r \ne 0\) ,而不是\(r>0.6\).
事实上，如果你希望推断"存在显著的强相关性"， 你的假设和备择假设应该是：
\(H_0: r\le 0.6\)
\(H_1: r> 0.6\)
而我们绝大多数统计分析软件，在计算相关系数的时候提供的p值一句的显著性检验的零假设都是\(H_0=0\), 只是对存在性进行了检验，并没有检验相关系数的强弱。
2.4 r=0.8, p<0.05 能说明存在显著的强相关吗？ No
不少人，论文甚至对统计分析很熟悉的人，也常见有人把一个很大的r值和一个显著的p值解释成存在显著的强相关性。这种错误解释的原因是对相关系数的显著性检验p值理解上的错误。 根据上一节的描述，p<0.05并不是对相关性强弱的显著性，而是相关性存在的显著性而已。如果从\(r=0.8, p<0.05\) 推断是否存在显著性的强相关性，可以用两种方法：
2.4.1 方法1：零假设显著性检验
用\(H_0:r \le 0.6\) , \(H_1:r>0.6\) 进行检验， 如果得到p<0.05 , 则拒绝零假设，推断存在显著的强相关性；否则，不能拒绝零假设 （！注意： 不能推断零假设成立）。
2.4.2 方法2：置信区间
因为样本的Pearson 相关系数并不服从正态分布，所以估计其置信区间，需要先进行Fisher -Z变换，将其样本分布正态化，然后根据z-分布计算出r在z-域的置信区间，再进行逆Z-变换， 得到r值的置信区间后，通过置信区间的范围来评估相关系数的强弱的显著性。我们以\(r=0.8, p=0.03, n=10\) 为例， 计算r的置信区间， R代码如下：
Code
library(DescTools)
#| echo: TRUE
pearsonr.ci<-function(r,ci,n){
alpha=1-ci
za=qnorm(alpha/2,0,1)
zr=0.5*log((1+r)/(1-r))
zL=zr-za*sqrt(1/(n-3))
zU=zr+za*sqrt(1/(n-3))
rL=(exp(2*zL)-1)/(exp(2*zL)+1)
rU=(exp(2*zU)-1)/(exp(2*zU)+1)
moeLU<-list("rL"=rL, "rU"=rU)
return(moeLU)
}
moe<-pearsonr.ci(0.8,0.95,10)
moe$rL
[1] 0.9507384
$rU
[1] 0.3432884可以看出r值的95%置信区间为\([0.34,0.95]\) , 这个范围与通常的弱相关（\(r<0.4\) ）, 中等强度线性相关（\(0.4<r<0.6\) ）, 强/极强相关（\(r>0.6\)）都存在完全或部分程度的重合，因此在样本量只有\(n=10\)的情况下，不能推断存在显著性的强相关。
类似的， 我们根据\(r=0.3, p=0.02\) 也不能推断存在显著的弱相关性，仍然要根据相关系数，样本量通过计算置信区间来推断。
3 Pearson相关系数易受异常值（outlier）影响
统计学里著名的数据集Anscombe， 展示了四组完全不同的数据，却具备相同的均值，中位数，标准差，相关系数。 这个数据集最常见的被用来演示在进行数据统计分析之前，进行数据可视化的重要性。我们以这个数据集的第四组数据为例，说明一个异常值会导致相关系数的很大变化， 比如anscombe数据集的相关系数达到\(r=0.82\), 其95%置信区间\([0.43, 0.9]\) , 甚至可以推断具有显著的中等以上的相关性强度了。
Code
library(datasets)
library(ggplot2)
datasets::anscombe x1 x2 x3 x4 y1 y2 y3 y4
1 10 10 10 8 8.04 9.14 7.46 6.58
2 8 8 8 8 6.95 8.14 6.77 5.76
3 13 13 13 8 7.58 8.74 12.74 7.71
4 9 9 9 8 8.81 8.77 7.11 8.84
5 11 11 11 8 8.33 9.26 7.81 8.47
6 14 14 14 8 9.96 8.10 8.84 7.04
7 6 6 6 8 7.24 6.13 6.08 5.25
8 4 4 4 19 4.26 3.10 5.39 12.50
9 12 12 12 8 10.84 9.13 8.15 5.56
10 7 7 7 8 4.82 7.26 6.42 7.91
11 5 5 5 8 5.68 4.74 5.73 6.89Code
lm4 <- lm(y4 ~ x4, data = anscombe)
cf<-coef(lm4)
eq = paste0("Dataset 4: y4 = ", round(cf[2],1), "*x4 +", round(cf[1],1))
pcc=paste0("Pearson r = ",round(cor(anscombe$x4,anscombe$y4),2),"\n p=",format(cor.test(anscombe$x4,anscombe$y4)$p.value,scientific = TRUE,digits=3))
p4 <- ggplot(anscombe) +
geom_point(aes(x4, y4), color = "darkorange", size = 1.5) +
scale_x_continuous(breaks = seq(0,20,2)) +
scale_y_continuous(breaks = seq(0,12,2)) +
expand_limits(x = 0, y = 0) +
labs(x = "x4", y = "y4",
title = eq ) +
theme_bw()+
annotate("text", label = pcc, x = 6, y = 10, size = 5, colour = "blue")
# geom_text()
#p4<-p4+
# geom_text()+
#
p4_fitted <- p4 + geom_abline(intercept = 3.0017, slope = 0.499, color = "blue")
p4_fitted
对Anscombe数据集感兴趣的，或者想了解这个数据集的其它几组数据的，可以参考Anscombe
鉴于相关系数很容易受到异常值的影响，我们一般建议在数据具有异常值（将数据可视化可以帮助发现），或者方差不具备齐性，即在不同自变量x值下，因变量y具有显著不同的方差特性，不建议使用Pearson 相关系数， 可以转而用Spearman’s \(\rho\), 或Kendall's \(\tau\), 比如下图展示的数据，就不适合用Pearson 相关系数进行分析：
Code
set.seed(10)
x=seq(from=1,to=50,by=1)
y=2*x+x*0.1*rnorm(50,0,10)
# plotting the regression line
#mod <- lm(y ~ x)
#cf <- coef(mod)
plot(x,y)
对于这样的数据，建议下面的方法之一：
将数据进行变换，基本满足方差齐性后，再计算Pearson相关系数。但是数据变换后，其结果的解释性将变差。
采用Spearman,或Kendall相关系数计算。