斯皮尔曼相关系数#
斯皮尔曼等级相关系数是一种衡量两个数据集之间单调关系强度的非参数度量。
考虑以下来自 [1] 的数据，该研究调查了不健康人肝脏中游离脯氨酸（一种氨基酸）和总胶原蛋白（一种常在结缔组织中发现的蛋白质）之间的关系。
下面的 x 和 y 数组记录了这两种化合物的测量值。这些观测值是成对的：每个游离脯氨酸测量值都与相同索引处的总胶原蛋白测量值取自同一个肝脏。
import numpy as np
# total collagen (mg/g dry weight of liver)
x = np.array([7.1, 7.1, 7.2, 8.3, 9.4, 10.5, 11.4])
# free proline (μ mole/g dry weight of liver)
y = np.array([2.8, 2.9, 2.8, 2.6, 3.5, 4.6, 5.0])
这些数据在 [2] 中使用斯皮尔曼相关系数进行了分析，这是一种对样本之间单调相关性敏感的统计量，在 scipy.stats.spearmanr 中实现。
from scipy import stats
res = stats.spearmanr(x, y)
res.statistic
np.float64(0.7000000000000001)
对于具有强正序数相关性的样本，该统计量的值趋于高（接近 1）；对于具有强负序数相关性的样本，其值趋于低（接近 -1）；对于具有弱序数相关性的样本，其值趋于小（接近零）。
该检验通过将统计量的观测值与零假设分布进行比较来执行：零假设分布是在总胶原蛋白和游离脯氨酸测量值独立这一零假设下推导出的统计量值分布。
对于此检验，可以将统计量进行转换，使得大样本的零假设分布是自由度为 len(x) - 2 的学生 t 分布。
import matplotlib.pyplot as plt
dof = len(x)-2 # len(x) == len(y)
dist = stats.t(df=dof)
t_vals = np.linspace(-5, 5, 100)
pdf = dist.pdf(t_vals)
fig, ax = plt.subplots(figsize=(8, 5))
def plot(ax): # we'll reuse this
ax.plot(t_vals, pdf)
ax.set_title("Spearman's Rho Test Null Distribution")
ax.set_xlabel("statistic")
ax.set_ylabel("probability density")
plot(ax)
plt.show()
该比较通过 p 值量化：p 值是零假设分布中与统计量观测值一样极端或更极端的值的比例。在双边检验中，如果统计量为正，则零假设分布中大于转换后统计量的元素以及小于观测统计量负值的元素都被视为“更极端”。
fig, ax = plt.subplots(figsize=(8, 5))
plot(ax)
rs = res.statistic # original statistic
transformed = rs * np.sqrt(dof / ((rs+1.0)*(1.0-rs)))
pvalue = dist.cdf(-transformed) + dist.sf(transformed)
annotation = (f'p-value={pvalue:.4f}\n(shaded area)')
props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
_ = ax.annotate(annotation, (2.7, 0.025), (3, 0.03), arrowprops=props)
i = t_vals >= transformed
ax.fill_between(t_vals[i], y1=0, y2=pdf[i], color='C0')
i = t_vals <= -transformed
ax.fill_between(t_vals[i], y1=0, y2=pdf[i], color='C0')
ax.set_xlim(-5, 5)
ax.set_ylim(0, 0.1)
plt.show()
res.pvalue
np.float64(0.07991669030889909)
如果 p 值“很小”——也就是说，从独立分布中抽样数据产生如此极端的统计量值的概率很低——这可以作为反对零假设并支持备择假设的证据：总胶原蛋白和游离脯氨酸的分布是不独立的。请注意，
反之则不成立；也就是说，该检验不用于提供支持零假设的证据。
被视为“小”值的阈值是在分析数据之前应做出的选择 [3]，同时应考虑假阳性（错误地拒绝零假设）和假阴性（未能拒绝一个错误的零假设）的风险。
小 p 值并非是大效应的证据；相反，它们只能提供“显著”效应的证据，这意味着它们在零假设下不太可能发生。
假设在进行实验之前，作者有理由预测总胶原蛋白和游离脯氨酸测量值之间存在正相关，并且他们选择根据单边备择假设来评估零假设的合理性：游离脯氨酸与总胶原蛋白存在正序数相关。在这种情况下，零假设分布中只有那些大于或等于观测统计量的值才被视为更极端。
res = stats.spearmanr(x, y, alternative='greater')
res.statistic
np.float64(0.7000000000000001)
fig, ax = plt.subplots(figsize=(8, 5))
plot(ax)
pvalue = dist.sf(transformed)
annotation = (f'p-value={pvalue:.6f}\n(shaded area)')
props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
_ = ax.annotate(annotation, (3, 0.018), (3.5, 0.03), arrowprops=props)
i = t_vals >= transformed
ax.fill_between(t_vals[i], y1=0, y2=pdf[i], color='C0')
ax.set_xlim(1, 5)
ax.set_ylim(0, 0.1)
plt.show()
res.pvalue
np.float64(0.03995834515444954)
请注意，t 分布提供了零假设分布的渐近近似；它只对有大量观测值的样本准确。对于小样本，执行置换检验可能更合适 [4]：在总胶原蛋白和游离脯氨酸独立的零假设下，每个游离脯氨酸测量值与任何总胶原蛋白测量值一起被观测到的可能性是均等的。因此，我们可以通过计算 x 和 y 之间所有可能的元素配对下的统计量来形成一个精确的零假设分布。
def statistic(x): # explore all possible pairings by permuting `x`
rs = stats.spearmanr(x, y).statistic # ignore pvalue
transformed = rs * np.sqrt(dof / ((rs+1.0)*(1.0-rs)))
return transformed
ref = stats.permutation_test((x,), statistic, alternative='greater',
permutation_type='pairings')
fig, ax = plt.subplots(figsize=(8, 5))
plot(ax)
ax.hist(ref.null_distribution, np.linspace(-5, 5, 26),
density=True)
ax.legend(['asymptotic approximation\n(many observations)',
f'exact \n({len(ref.null_distribution)} permutations)'])
plt.show()
ref.pvalue
np.float64(0.04563492063492063)