如何计算 Pearson 相关矩阵并仅保留重要值？

Posted 2023-02-16

【中文标题】如何计算 Pearson 相关矩阵并仅保留重要值？

【英文标题】：How can I compute the Pearson correlation matrix and retain only significant values?

【发布时间】：2020-04-12 10:39:57

【问题描述】：

我有一个 4×3 矩阵 X，并希望通过计算 X 的所有 3 个可能列组合之间的相关性来形成 3×3 Pearson 相关矩阵 C。但是，C 中对应于统计上不显着的相关性的条目应设置为零。

我知道如何在scipy.stats 中使用pearsonr 获得成对相关性和显着性值。例如，

import numpy as np
from scipy.stats.stats import pearsonr

X = np.array([[1, 1, -2], [0, 0, 0], [0, .2, 1], [5, 3, 4]])
pearsonr(X[:, 0], X[:, 1])

返回(0.9915008164289165, 0.00849918357108348)，X 的第一列和第二列之间的相关性约为 0.9915，p 值为 0.0085。

我可以使用嵌套循环轻松获得所需的矩阵：

将 C 预填充为 3×3 零矩阵。 嵌套循环的每一遍将对应于X 的两列。如果 p 值小于或等于我的阈值，例如 0.01，则对应于这对列的 C 条目将设置为成对相关。

我想知道是否有更简单的方法。我知道在 Pandas 中，我可以在基本上一行中创建相关矩阵 C：

import pandas as pd

df = pd.DataFrame(data=X)
C_frame = df.corr(method='pearson') 
C = C_frame.to_numpy()

有没有办法在没有循环的情况下获取 p 值的矩阵或数据框P？如果是这样，如果P 中的相应p 值超过我的阈值，我如何将C 的每个条目设置为零？

【问题讨论】：

C_frame.where(C_frame>0.99)?

@QuangHoang。这根本不是一回事

***.com/questions/52741236/… 是相关的。高度投票的答案显示了如何使用 method 参数返回 p-values 而不是相关系数。您可以将其用于mask 您的df.corr() 结果。虽然它仍然是一个循环......

仅供参考：***.com/a/38946645/1217358

【参考方案1】：

查看pearsonr 的文档会发现用于计算相关性的公式。使用向量化得到矩阵每一列之间的相关性应该不会太难。

虽然您可以使用 pandas 计算 C 的值，但我将展示整个过程的纯 numpyan 实现。

首先，计算 r 值：

X = np.array([[1,  1, -2],
              [0,  0,  0],
              [0, .2,  1],
              [5,  3,  4]])
n = X.shape[0]

X -= X.mean(axis=0)
s = (X**2).sum(axis=0)
r = (X[..., None] * X[..., None, :]).sum(axis=0) / np.sqrt(s[:, None] * s[None, :])

考虑到 scipy 中存在 beta 分布，计算 p 值变得很简单。直接取自文档：

dist = scipy.stats.beta(n/2 - 1, n/2 - 1, loc=-1, scale=2)
p = 2 * dist.cdf(-abs(r))

您可以使用您的阈值从p 轻松制作一个蒙版，并将其应用于r 以制作C：

mask = (p <= 0.01)
C = np.zeros_like(r)
C[mask] = r[mask]

更好的选择可能是就地修改您的r：

r[p > 0.1] = 0

函数形式：

def non_trivial_correlation(X, threshold=0.1):
    n = X.shape[0]
    X = X - X.mean(axis=0) # Don't modify the original
    x = (X**2).sum(axis=0)
    r = (X[..., None] * X[..., None, :]).sum(axis=0) / np.sqrt(s[:, None] * s[None, :])
    p = 2 * scipy.stats.beta(n/2 - 1, n/2 - 1, loc=-1, scale=2).cdf(-abs(r))
    r[p > threshold] = 0
    return r