4. 走进 Combat 的数学内核：位置、尺度与贝叶斯魔法

好，我们开始啃硬骨头了。ComBat 的核心，说白了就三件事：找到批次间的差异、量化这个差异、然后把它消除掉。数学上，它把批次效应拆成了两个部分：一个是“整体偏移”，就是位置参数；另一个是“波动幅度”，就是尺度参数。

我个人习惯把这套流程想象成调音师的工作。一台钢琴（你的数据），在不同音乐厅（不同批次）里，音准和共鸣都不一样。调音师得先听出每个厅的“偏移”和“混响”，再统一调回标准。ComBat 做的就是这个活儿。

4.1 位置参数（μ）与尺度参数（δ）：批次效应的两个维度

在 ComBat 的数学世界里，我们把每个特征（比如一个蛋白质的表达强度）拆开来建模。它的经典模型长这样：

Y_{ijg} = α_g + Xβ_g + γ_{ig} + δ_{ig}ε_{ijg}

我来一个个拆解给你看：

  • Y_{ijg}: 第 i 个批次、第 j 个样本、第 g 个特征的原始表达值。这就是你拿到的原始数据。
  • α_g: 特征 g 的全局基线水平。相当于钢琴的中音 C 的基准频率。
  • Xβ_g: 生物学协变量的影响矩阵。比如疾病组 vs 对照组，这部分差异是我们要保留的，不能给“校掉”了。
  • γ_{ig}: 批次 i 对特征 g 的 位置偏移 (location shift)。这就是我前面说的“音准偏差”。
  • δ_{ig}: 批次 i 对特征 g 的 尺度缩放 (scale adjustment)。这是“混响大小”。
  • ε_{ijg}: 随机误差项。通常假设它服从标准正态分布 N(0, σ_g²)。

你看，γ 负责把数据整体拉高或拉低，δ 负责把数据的离散程度放大或缩小。校准的目标就是：估计出 γ 和 δ，然后把它们从原始数据中“剥离”出去。

核心校准公式（我习惯称它为“去批次”操作）：
校正后的数据 Z_{ijg} = (Y_{ijg} - α_g - Xβ_g - γ*_{ig}) / δ*_{ig}
这里的 γ* 和 δ* 就是我们估计出来的批次效应参数。除一下 δ，别忘了。

4.2 为什么用经验贝叶斯？

你可能要问了：“这么简单，直接拿每个批次的均值和方差算 γ 和 δ 不就行了？” 嗯，我一开始也是这么想的。

但现实很骨感。质谱数据有个特点：特征多（成千上万个蛋白/肽段），批次少（通常只有两三个批次）。你直接用每个批次的数据去估算某个特征的 δ，会发现方差极大，估计极其不稳定。尤其是一些表达量低的特征，信噪比本来就差，估计出的 δ 可能直接炸掉，导致校正后数据出现极端值。

这就是 ComBat 的高明之处——它引入了经验贝叶斯 (Empirical Bayes, EB) 框架。

说白了，EB 就是“借用力量”（borrowing strength）。在估算某一个特征的 γ 和 δ 时，它会偷看其他成千上万个特征的信息，利用全局规律来约束单个特征的估计值，让结果更稳健。我遇到过一些项目，不用 EB 直接做位置-尺度校正，结果校正后的 PCA 图上，批次效应反而变得更乱了。这坑我踩过，所以印象很深。

经验贝叶斯的流程可以抽象为下图，我习惯把它看作一个“闭环反馈”系统：

经验贝叶斯框架下的ComBat参数估计流程 输入标准化数据 Y_{ijg} - α_g - Xβ_g 标准矩估计 γ̂_ig, δ̂_ig (粗糙) 先验分布的设定 γ ~ N(μγ, σ²γ) δ² ~ InvGamma(λ, θ) 计算后验分布（EB） 收缩估计 γ*_ig, δ*_ig 输出校正数据 Z_{ijg} = (Y_ijg - γ*) / δ* 超参数拟合 经验贝叶斯的核心贡献： 用单个特征的矩估计结果去拟合全局先验分布的“超参数”，再用这个先验去“拉回”那些不可靠的单个估计值。
个人经验： 我在处理 DIA 质谱数据时，发现一些低丰度蛋白的 CV 值（变异系数）高达 80%。直接做位置-尺度校正，这些蛋白的 δ 估计值奇高。使用 EB 后，它们的 δ* 会被大幅度“收缩”到全局中位数附近，校正后的CV值稳定在 20% 以内。

4.3 先验分布怎么设定？——数学家的“押注”

好，现在到了 EB 框架最关键的一步：先验分布。ComBat 的作者为 γ 和 δ² 选择了 共轭先验。你可能会问，为什么选共轭的？说白了，就是方便计算后验分布，能用数学公式直接算出来，不用跑 MCMC 模拟那么费劲。

具体设定如下：

参数 先验分布 超参数 直观理解
位置参数 γ 正态分布 N(μ_γ, σ²_γ) μ_γ, σ²_γ 假设大多数批次间的整体偏移集中在某个均值附近，不会偏离太离谱。
尺度参数 δ² 逆伽马分布 InvGamma(λ, θ) λ (形状), θ (尺度) 假设大多数特征的方差缩放因子 δ² 服从一个偏态分布，不会太大也不会太小。

这里有两个细节值得注意：

  1. 超参数的估计：μ_γ、σ²_γ、λ、θ 这些被叫做“超参数”，它们本身也是从数据里估算出来的。怎么估？就是用我们前面拿到的那一堆粗糙的矩估计值 (γ̂, δ̂) 去拟合。这就是“经验”贝叶斯这个名字的由来——先验不是拍脑门定的，是从经验数据里学出来的。
  2. 为什么用逆伽马：δ² 是方差，它只能是正数。正态分布不适合。逆伽马分布是定义在正数上的，而且它是伽马分布共轭于正态分布方差参数的先验。选它，数学上最干净，解析解最优美。嗯，这里要注意，ComBat 论文里假设的是 δ² 服从逆伽马，而不是 δ。我早期看代码时在这里卡了很久。
避坑指南： 我曾经遇到过一批质谱数据，其中有一部分特征的信号完全缺失（全是0）。如果用这些特征去拟合逆伽马先验，会严重拉偏分布的形状。我的经验是：在运行ComBat前，务必先过滤掉那些缺失率超过 50% 的特征，或者用 KNN 做一下填补。否则，先验分布会学到一个非常扭曲的形状，导致校正效果变差。

有了这些先验，ComBat 就能计算出每个特征 γ 和 δ² 的后验分布。然后取后验分布的期望或中位数，作为最终的 γ* 和 δ*。这个“收缩”过程，就是经验贝叶斯的核心魔法——它用全局信息支撑了个体估计，解决了批次少、特征多的统计难题。

← 上一章 📖 返回目录 下一章 →