在python中获取相关系数的统计差异

Question

为了获得python中两个数组之间的相关性，我正在使用：

from scipy.stats import pearsonr
x, y = [1,2,3], [1,5,7]
cor, p = pearsonr(x, y)

但是，如文档中所述，从

pearsonr()

返回的p值仅对大于500的数据集有意义。那么如何获得对于小数据集来说合理的p值？

我的临时解决方案：

在阅读了线性回归之后，我想出了自己的小脚本，它基本上使用 Fischer 变换 来获取 z 分数，从中计算 p 值：

import numpy as np
from scipy.stats import zprob
n = len(x)
z = np.log((1+cor)/(1-cor))*0.5*np.sqrt(n-3))
p = zprob(-z)

它有效。但是，我不确定

pearsonr()

给出的 p 值是否更合理。是否有一个Python模块已经具有此功能？我在 SciPy 或 Statsmodels 中找不到它。

编辑以澄清：

我的示例中的数据集经过简化。我的真实数据集是两个包含 10-50 个值的数组。

Answer 1

scipy.stats.pearsonr

现在接受

scipy.stats.PermutationMethod

的实例作为其

method

参数传递。这使您可以对原假设进行精确检验，即

x

和

y

中的观测值是从独立（因此不相关）分布

X

和

Y

中得出的。在零假设下，

y

中元素的所有可能排列都是同样可能的。通过计算所有可能排列下的统计量，我们可以计算获得与实际观察到的极端统计量一样极端的检验统计量的确切概率。

对于您的数据：

from scipy import stats
x, y = [1, 2, 3], [1, 5, 7]
pearsonr(x, y, method=stats.PermutationMethod())
# PearsonRResult(statistic=0.9819805060619655, pvalue=0.3333333333333333)

对于较大的样本，可能的排列数量太大，因此进行随机检验。

有关 p 值仅对如此大的样本有意义的注释已在最新版本的文档中删除。对于大小为 500 的样本，p 值非常匹配（取决于随机检验的随机性质），尽管从技术上讲，两个版本的检验的具体原假设略有不同。

import numpy as np
from scipy import stats
rng = np.random.default_rng(2348923589435)
x, y = rng.random((2, 20))
stats.pearsonr(x, y)
# PearsonRResult(statistic=0.17895011346631318, pvalue=0.45031837813125064)
stats.pearsonr(x, y, method=stats.PermutationMethod(random_state=rng))
# PearsonRResult(statistic=0.17895011346631318, pvalue=0.4412)

然而，不需要如此大的样本就能获得类似的一致性。例如，对于每个样本大小 4-30，我们可以生成由

pearsonr

（默认

method

）生成的 p 值为 0.05 的数据，然后与排列测试生成的 p 值进行比较。

# it would be better to compare the null distributions directly,
# but this is easier without digging into the implementations
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats, optimize

rng = np.random.default_rng(2348923589435)
ns = np.arange(4, 31)
diffs = []
for n in ns:
    x = rng.random(size=n)
    d = rng.normal(size=n)
    res = optimize.root_scalar(lambda c: stats.pearsonr(x, x + c*d).pvalue - 0.05, 
                               bracket=(0, 10))
    y = x + res.root*d
    _, p1 = stats.pearsonr(x, y)
    _, p2 = stats.pearsonr(x, y, method=stats.PermutationMethod())
    diffs.append(np.abs(p2 - p1))

plt.plot(ns, diffs)

对于任何大于 10 的样本，两者似乎非常一致。 有关更多信息，请参阅重采样和蒙特卡罗方法教程。