第3章 批次效应的统计检验：ANOVA、Kruskal-Wallis与置换检验

3.1 为什么要先做统计检验？

你想想看，做完PCA或者t-SNE图，看到不同批次在图上分开了。那是不是一定就有批次效应？

不一定。视觉这个东西容易骗人。我个人习惯，不管图上看起来怎么样，都会跑一遍统计检验来核实一下。

我在项目中遇到过好几次，PCA图上样本混在一起，看着没问题。结果ANOVA一跑，p值显著得很。这说明什么？说明批次效应藏在高维空间里，降维给漏掉了。

反过来也有，图上分开得很明显，但检验下来根本不显著。可能是样本量太小，或者某个离群点把图给拉偏了。

所以，我的标准流程里，统计检验是必须的。根据数据情况，从ANOVA、Kruskal-Wallis，还有置换检验里挑一个。

3.2 ANOVA（方差分析）——参数检验利器

ANOVA检验的是多个批次的均值是否相等。它要求数据近似正态分布，并且组间方差差不多。

在组学数据里，我们通常对每个基因（或代谢物、蛋白）独立做ANOVA。

代码其实很简单，用scipy.stats就能搞定：

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

# 假设data是数据框，行为特征，列为样本
# batch_labels是批次标签列表
# 对第一个特征做ANOVA
feature = data.iloc[0, :]
groups = [feature[batch_labels == b] for b in np.unique(batch_labels)]
f_stat, p_value = f_oneway(*groups)

print(f"F统计量: {f_stat}")
print(f"P值: {p_value}")

跑完所有特征后，你会得到一大堆p值。这时候要注意——千万别忘了多重假设检验校正。

import statsmodels.stats.multitest as smm

# 假设p_values是所有特征的原始p值列表
rejected, corrected_p_values = smm.multipletests(p_values, method='fdr_bh')

3.3 Kruskal-Wallis检验——非参数替代方案

组学数据有个特点——经常不符合正态分布。比如代谢组学，受生理波动影响大，离群值多。这时候ANOVA就不太合适了。

Kruskal-Wallis检验是ANOVA的非参数替代。它不关心原始值，只关心秩（rank），对异常值稳健得多。

代码几乎一样，换个函数的事：

from scipy.stats import kruskal

feature = data.iloc[0, :]
groups = [feature[batch_labels == b] for b in np.unique(batch_labels)]
h_stat, p_value = kruskal(*groups)

print(f"H统计量: {h_stat}")
print(f"P值: {p_value}")

我个人经验，做代谢组学数据分析时，我会优先用Kruskal-Wallis。代谢组数据经常有异常值，KW检验不容易被那几个离群点带偏。

嗯，为什么会这样？因为它是基于秩的。1万和1亿，在KW看来可能只是第1名和第2名的区别。不会因为数值大就被当成异常。

3.4 置换检验——万金油方法

如果你不确定数据符不符合参数检验的假设，或者样本量太小，常规检验不好使，那怎么办？

置换检验就是万金油。它不依赖任何理论分布，而是通过反复打乱批次标签，生成统计量的零分布。然后看你的观测值在零分布里处在什么位置。

说白了，就是反复问自己：“如果完全没有批次效应，随机打乱标签，结果比现在更极端的概率有多大？”

import numpy as np
from scipy.stats import f_oneway

def permutation_test(data, labels, n_permutations=5000):
    observed_f, _ = f_oneway(*[data[labels == b] for b in np.unique(labels)])
    count = 0
    all_labels = labels.copy()
    np.random.seed(42)
    for _ in range(n_permutations):
        np.random.shuffle(all_labels)
        perm_f, _ = f_oneway(*[data[all_labels == b] for b in np.unique(all_labels)])
        if perm_f >= observed_f:
            count += 1
    p_value = (count + 1) / (n_permutations + 1)  # +1避免p=0
    return observed_f, p_value

f_stat, p_val = permutation_test(data.iloc[0, :].values, batch_labels)
print(f"置换检验P值: {p_val}")

置换次数越多，p值精度越高。我一般至少跑5000次，重要分析上到10000次。代价就是计算量会大一些。

3.5 三种方法对比总结

为了让你看得更清楚，我列了个对比表：

下面这是我自己总结的决策流程，照着走一般不会错：

3.6 总结与实用建议