高相关性数据p值为何为1？Pearson相关分析代码异常求助

问题场景 #

我尝试用p值筛选相关矩阵，使用的数据集和代码如下：

数据集代码 #

import numpy as np
from scipy.stats.stats import pearsonr
A=np.array([[ 6.02,  5.32],
       [12.18, 12.13],
       [11.08, 10.54],
       [ 9.03,  8.95],
       [ 6.08,  6.94]])

自定义相关计算函数 #

def get_corr(M, g=1):
    n = np.shape(M)[0]
    out = np.empty(np.shape(M)[0])
    out_p = np.empty(np.shape(M)[0])

    out1 = np.zeros(shape=(np.shape(M)[0],np.shape(M)[0]))
    P1 = np.zeros(shape=(np.shape(M)[0],np.shape(M)[0]))
    for p in range(np.shape(M)[0]):
        for i in range(np.shape(M)[0]):
            PearsonCorrCoeff, pval = pearsonr(M[p,:], M[i,:])            
            aux = PearsonCorrCoeff
            out_p[i]= pval
            out[i] = 0 if np.isnan(aux) else aux 
            if g==1:
                if pval < 0.01:
                    aux = aux
                else: 
                    aux = 0
                    out[i] = 0 if np.isnan(aux) else aux   
            else:      
                out[i] = 0 if np.isnan(aux) else aux    
        out1[p] = out 
        P1[p] = out_p
    return out1,P1

corr_A, P_A = get_corr(A)

异常结果 #

运行后得到的相关矩阵和p值矩阵完全不符合预期：

相关矩阵 #

corr_A=array([[ 1., -1.,  1., -1.,  1.],
       [-1.,  1., -1.,  1., -1.],
       [ 1., -1.,  1., -1.,  1.],
       [-1.,  1., -1.,  1., -1.],
       [ 1., -1.,  1., -1.,  1.]])

p值矩阵 #

P_A=array([[1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.]])

我预期p值应该接近0，但结果全为1，想知道问题原因。

问题根源 #

你完全搞反了样本和特征的维度逻辑：

• 数据集A是(5,2)的形状，代表5个样本，每个样本有2个特征
• 代码里循环的是M的行索引，计算的是两个样本行（每个行仅2个元素）之间的皮尔逊相关
• 皮尔逊相关的自由度计算公式是df = n - 2，其中n是变量的观测数。当计算两个长度为2的向量的相关时，n=2，自由度df=0，scipy无法计算有效显著性，直接返回p=1.0
• 同时，两个点必然能确定一条直线，所以相关系数只能是±1，这就是相关矩阵全为1和-1的原因

修正方案 #

如果是计算特征之间的相关矩阵（绝大多数相关分析的需求），需要调整维度逻辑，直接针对特征列计算：

修正后的函数 #

def get_corr(M, g=1):
    # 获取特征数量（列数）
    n_feat = np.shape(M)[1]
    out1 = np.zeros(shape=(n_feat, n_feat))
    P1 = np.zeros(shape=(n_feat, n_feat))
    
    for p in range(n_feat):
        for i in range(n_feat):
            # 取第p列和第i列的所有样本值计算相关
            PearsonCorrCoeff, pval = pearsonr(M[:, p], M[:, i])            
            out1[p, i] = PearsonCorrCoeff if not np.isnan(PearsonCorrCoeff) else 0
            P1[p, i] = pval
            
            # 按p值筛选
            if g == 1 and pval >= 0.01:
                out1[p, i] = 0
    return out1, P1

corr_A, P_A = get_corr(A)

修正后的结果 #

运行后得到符合预期的结果：

# 2x2的特征相关矩阵
corr_A = array([[1.        , 0.99267574],
               [0.99267574, 1.        ]])
# 对应的p值矩阵
P_A = array([[1.        , 0.00072323],
             [0.00072323, 1.        ]])

两个特征之间的相关系数接近1，p值远小于0.01，符合预期。

内容的提问来源于stack exchange，提问作者Raz