nan-safe-correlation

NaN 安全的相关性计算

概述

存在缺失值时，跨多个特征（基因、蛋白质、变体）计算相关性容易出错。最常见的错误是使用批量矩阵快捷方法，这类方法会悄无声息地错误处理 NaN，导致相关性数值不准确。本指南介绍正确的逐特征配对计算、退化输入过滤和性能优化。

关键概念

配对删除与列表删除

• 配对删除：对每一对特征，仅移除任一数值为 NaN 的样本。每个特征对都能使用最多的可用数据。
• 列表删除：移除任一特征中存在 NaN 的所有样本。这会浪费有效数据；如果缺失并非完全随机，还会造成结果偏差。
• 规则：逐特征计算相关性时，始终使用配对删除。

批量矩阵快捷方法为何会失效

不同特征在各样本中的缺失值模式各不相同。批量方法对这些情况的处理并不一致：

方法	问题
DataFrame.rank() 然后使用 corrwith()	rank() 会为 NaN 分配秩；corrwith() 可能会全局或按列不一致地删除数据
DataFrame.corrwith(method='spearman')	实现因 pandas 版本而异；可能使用列表删除
对秩数据使用 np.corrcoef	只要有一个值缺失，就会将 NaN 传播到整个结果

错误计算的影响

• 相关性可能偏差 0.01–0.05 或更多
• 接近阈值（例如 0.6）的特征可能被错误分类
• 每个特征的有效样本数未知（可能在不知情的情况下使用少于预期的样本）

退化输入

会产生未定义或不稳定相关性的特征：

类型	描述	影响
常量特征	所有数值都相同（方差 = 0）	相关性未定义（除以零）
近常量特征	方差非常低	相关性在数值上不稳定
有效值过少	移除 NaN 后，有效配对少于 min_valid	统计结果不可靠
过滤后仅剩单个值	移除 NaN 后只剩一个唯一值	相关性未定义

决策框架

Do you have missing values (NaN) in your feature matrix?
├── No NaN at all → Bulk methods are safe (corrwith, np.corrcoef)
└── Yes, NaN present
    ├── Same NaN pattern across all features? → Listwise deletion is acceptable
    └── Different NaN patterns per feature (typical)
        ├── < 10,000 features → Per-feature loop with scipy.stats.spearmanr
        └── > 10,000 features → Parallelized per-feature loop (joblib)

场景	推荐方法	理由
无缺失数据	DataFrame.corrwith()	速度快；无 NaN 时结果正确
NaN 稀疏，特征数 < 10K	逐特征 spearmanr 循环	正确执行配对删除，速度可接受
NaN 稀疏，特征数 > 10K	并行化的逐特征循环	正确性相同，可利用多核扩展
NaN 密集（缺失率 > 50%）	逐特征循环 + 严格的 min_valid	许多特征会被跳过；报告跳过数量
NaN 模式一致	列表删除 + 批量方法	如果所有特征共享相同的 NaN 行，则配对删除 = 列表删除

最佳实践

1. 始终在分析前打印 NaN 摘要：报告 NaN 总数、包含 NaN 的特征，以及每个特征的 NaN 分布。这有助于记录数据质量，并提醒你注意严重的缺失模式。

2. 在循环中对每个特征使用 scipy.stats.spearmanr：这是唯一能确保对每个特征分别正确删除成对 NaN 的方法。

3. 设置有效样本对数下限（min_valid）：默认值设为 10。删除 NaN 后有效样本对数不足的特征会产生不可靠的相关性，应以 NaN 跳过。

4. 计算相关性前过滤退化输入：在相关性循环前移除常量特征、近似常量特征和 NaN 过多的特征。这样可以避免未定义结果并加快计算。

5. 在输出中记录每个特征的 n_valid：每个特征的有效样本对数各不相同。应将其与 rho 和 p 值一并报告，以便后续分析评估结果的可靠性。

6. 报告跳过或过滤的特征数量：特征无声丢失是造成困惑的常见原因。始终打印过滤掉的退化特征数和因数据不足而跳过的特征数。

7. 对大型数据集使用并行处理：对于超过 10,000 个特征，使用 joblib 将逐特征循环分配到多个核心。逐特征计算非常适合并行处理。

常见陷阱

1. 存在 NaN 时先批量排序再计算相关性：先调用 df.rank()，再调用 corrwith()，会在不提示的情况下错误处理 NaN，导致相关性计算不正确。

   • 避免方法：存在 NaN 时，始终对每个特征使用 scipy.stats.spearmanr。

2. 假设所有特征的样本数相同：不同特征的 NaN 模式各异，因此每次相关性计算使用的样本数也不同。

   • 避免方法：跟踪并报告每个特征的 n_valid。

3. 不对退化输入进行过滤：常量或近似常量特征会产生未定义的相关性或除以零警告，可能在不提示的情况下导致结果出错。

   • 避免方法：在相关性循环前运行 filter_degenerate()。

4. NaN 模式不同时使用按行删除：按行删除会移除任何特征中含有 NaN 的行，可能导致丢弃大部分数据。

   • 避免方法：使用成对删除（逐特征移除 NaN）。

5. 忽略 NaN 摘要步骤：跳过数据质量报告后，就无法判断 NaN 模式是否严重到足以影响结果。

   • 避免方法：始终在计算相关性前打印 NaN 摘要。

6. 将 min_valid 设得过低：有效样本对数少于约 10 时，Spearman 相关性不可靠，p 值也没有意义。

   • 避免方法：使用 min_valid >= 10；对于高维研究，应提高该值。

工作流

1. 步骤 1：打印 NaN 摘要

   • 报告数据集形状、NaN 总数、包含 NaN 的特征数，以及每个特征的平均和最大 NaN 数量
   • 决策点：如果总体 NaN 超过 50%，应在继续前重新评估数据质量

2. 步骤 2：过滤退化特征

   • 移除常量特征（nunique < 3）
   • 移除 NaN 过多的特征（非 NaN 值占比 < 50%）
   • 报告移除的特征数

3. 步骤 3：计算每个特征的相关性

   • 使用 scipy.stats.spearmanr 遍历各个特征
   • 对每个特征分别执行成对 NaN 删除
   • 跳过有效样本对数少于 min_valid 的特征（记录为 NaN）

4. 步骤 4：整理并报告结果

   • 创建包含每个特征的 rho、p 值和 n_valid 的 DataFrame
   • 报告跳过的特征数
   • 确认没有数据无声丢失（输入特征数 = 输出特征数 + 跳过数 + 过滤数）

参考实现

from scipy.stats import spearmanr
import numpy as np
import pandas as pd

def nan_summary(df):
    """Print NaN summary before correlation analysis."""
    print(f"Dataset shape: {df.shape}")
    print(f"Total NaN: {df.isna().sum().sum()}")
    print(f"Features with any NaN: {(df.isna().any()).sum()}")
    print(f"NaN per feature (mean): {df.isna().sum().mean():.1f}")
    print(f"NaN per feature (max): {df.isna().sum().max()}")

def filter_degenerate(df, min_unique=3, min_nonnan_frac=0.5):
    """Remove degenerate features before correlation analysis.

    Args:
        df: DataFrame (samples x features)
        min_unique: Minimum number of unique non-NaN values required
        min_nonnan_frac: Minimum fraction of non-NaN values required

    Returns:
        Filtered DataFrame, count of removed features
    """
    n_samples = len(df)
    keep = []
    for col in df.columns:
        values = df[col].dropna()
        if len(values) < n_samples * min_nonnan_frac:
            continue
        if values.nunique() < min_unique:
            continue
        keep.append(col)
    removed = len(df.columns) - len(keep)
    print(f"Filtered {removed} degenerate features out of {len(df.columns)}")
    return df[keep], removed

def pairwise_spearman(df_x, df_y, min_valid=10):
    """Compute per-feature Spearman correlation with pairwise NaN removal.

    Args:
        df_x: DataFrame (samples x features), aligned with df_y
        df_y: DataFrame (samples x features), same shape as df_x
        min_valid: Minimum number of valid (non-NaN) pairs required

    Returns:
        DataFrame with columns: rho, pvalue, n_valid
    """
    nan_summary(df_x)
    nan_summary(df_y)

    results = []
    for feature in df_x.columns:
        x = df_x[feature].values
        y = df_y[feature].values
        mask = ~(np.isnan(x) | np.isnan(y))
        n_valid = mask.sum()
        if n_valid < min_valid:
            results.append({'feature': feature, 'rho': np.nan,
                          'pvalue': np.nan, 'n_valid': n_valid})
            continue
        rho, pval = spearmanr(x[mask], y[mask])
        results.append({'feature': feature, 'rho': rho,
                       'pvalue': pval, 'n_valid': n_valid})

    result_df = pd.DataFrame(results).set_index('feature')
    skipped = result_df['rho'].isna().sum()
    if skipped > 0:
        print(f"Skipped {skipped} features with < {min_valid} valid pairs")
    return result_df

反模式

# WRONG: Bulk rank-then-correlate
ranked_x = df_x.rank()
ranked_y = df_y.rank()
corrs = ranked_x.corrwith(ranked_y)

# WRONG: Bulk corrwith with method parameter
corrs = df_x.corrwith(df_y, method='spearman')

# WRONG: numpy corrcoef on ranked arrays (propagates NaN)
corrs = np.corrcoef(df_x.rank().values.T, df_y.rank().values.T)

性能优化（超过 10,000 个特征）

from joblib import Parallel, delayed

def parallel_spearman(df_x, df_y, min_valid=10, n_jobs=4):
    """Parallelized per-feature Spearman correlation."""
    def compute_one(feature):
        x = df_x[feature].values
        y = df_y[feature].values
        mask = ~(np.isnan(x) | np.isnan(y))
        n = mask.sum()
        if n < min_valid:
            return feature, np.nan, np.nan, n
        rho, pval = spearmanr(x[mask], y[mask])
        return feature, rho, pval, n

    results = Parallel(n_jobs=n_jobs)(
        delayed(compute_one)(f) for f in df_x.columns
    )
    return pd.DataFrame(
        results, columns=['feature', 'rho', 'pvalue', 'n_valid']
    ).set_index('feature')

延伸阅读

• scipy.stats.spearmanr 文档 -- 官方 API 参考
• pandas corrwith 文档 -- 已知的 NaN 处理行为
• 成对删除与列表删除 -- 《缺失数据的灵活插补》中的章节

相关技能

• statistical-analysis -- 统计检验的一般选择方法与假设检查
• degenerate-input-filtering -- 在进行任何统计检验前筛除无信息数据的综合指南
• scikit-learn-machine-learning -- 特征选择与预处理流水线