nan-safe-correlation

NaN セーフな相関計算

概要

欠損値がある状態で多数の特徴量（遺伝子、タンパク質、変異）の相関を計算すると、エラーが発生しやすくなります。最もよくある間違いは、NaN を暗黙的に誤処理して不正な相関値を生成する一括行列処理のショートカットを使うことです。このガイドでは、特徴量ごとの正しいペアワイズ計算、退化した入力のフィルタリング、パフォーマンスの最適化について説明します。

主な概念

ペアワイズ削除とリストワイズ削除

• ペアワイズ削除：特徴量の各ペアについて、どちらかの値が NaN であるサンプルのみを除外します。各特徴量で利用可能なデータを最大限に使用します。
• リストワイズ削除：いずれかの特徴量に NaN があるサンプルをすべて除外します。有効なデータが無駄になり、欠損が完全にランダムでない場合は結果にバイアスが生じます。
• ルール：特徴量ごとの相関には、必ずペアワイズ削除を使用してください。

一括行列処理のショートカットが機能しない理由

サンプル間で、特徴量ごとに欠損値のパターンが異なります。一括処理のメソッドでは、この扱いに一貫性がありません。

メソッド	問題
DataFrame.rank() の後に corrwith()	rank() は NaN のランクを割り当てます。corrwith() は全体単位または列単位で一貫性なく除外する場合があります。
DataFrame.corrwith(method='spearman')	実装は pandas のバージョンによって異なり、リストワイズ削除を使用する場合があります。
ランク付けしたデータに対する np.corrcoef	値が 1 つでも欠損していると、結果全体に NaN が伝播します。

不正な計算による影響

• 相関値が 0.01～0.05、またはそれ以上変動することがあります。
• 閾値（例：0.6）付近の特徴量が誤分類されることがあります。
• 特徴量ごとの有効なサンプル数が分からず、想定より少ないサンプル数が暗黙的に使われる場合があります。

退化した入力

相関が未定義または不安定になる特徴量：

種類	説明	影響
定数特徴量	すべての値が同じ（分散 = 0）	相関は未定義（ゼロ除算）
ほぼ定数の特徴量	分散が非常に小さい	相関が数値的に不安定
有効値が少なすぎる	NaN の削除後に、有効なペアが min_valid 未満	統計的に信頼できない
フィルタリング後に値が 1 つ	NaN の削除後に、異なる値が 1 つしか残らない	相関は未定義

判断の指針

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 未満の場合、スピアマン相関の信頼性が低くなり、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 処理の動作
• ペアワイズ削除とリストワイズ削除 -- 『Flexible Imputation of Missing Data』の章

関連スキル

• statistical-analysis -- 統計検定の選択と前提条件の確認に関する一般的なガイド
• degenerate-input-filtering -- 統計検定を行う前に情報量の少ないデータをフィルタリングするための包括的なガイド
• scikit-learn-machine-learning -- 特徴選択と前処理パイプライン