第4节：SVA模型推导——可替代变量的两步估计法

好，咱们进入正题。SVA（Surrogate Variable Analysis）的核心思想，说白了就是：

数据里既藏着真正的生物学信号，也混着看不见的批次效应。我们没办法直接测量这些效应，那就造一些“可替代变量”来代表它们。

我个人习惯把SVA的流程拆成两步：
第一步：猜一下到底需要多少个可替代变量。
第二步：把这些变量具体算出来。

这两步的顺序不能乱。先确定数量，再算具体值。下面逐个拆解。

4.1 第一步：估计可替代变量的数量（K值）

SVA的第一步不走寻常路，它用的是平行分析（Parallel Analysis）。你可能在其他地方也见过这个名字，比如在因子分析里。

平行分析的基本思路是这样的：
把原始数据的特征值，和随机打乱后的数据特征值做个对比。

具体怎么操作？我整理成几个步骤：

1. 对原始数据做SVD（奇异值分解）。
   得到一组真实特征值，我们记为 λ_j（j = 1 到 m）。
2. 把数据矩阵的行（样本）随机打乱。
   这样做会把所有样本间的结构关系全部破坏掉，只留下独立同分布的噪声。
3. 对打乱后的数据再做一次SVD。
   得到一组“虚拟”特征值，记为 λ_j⁽⁰⁾。
4. 比较。
   如果某个 λ_j 明显大于 λ_j⁽⁰⁾，说明这个方向上有信号存在，不是纯噪声。

这里有公式：

K = max{ j | λ_j > quantile(λ_j^(0), α) }

其中 α 是显著性水平，通常取0.05。我一般会把这个过程重复 B 次（比如 B=100），然后取中位数或均值。这样估计出来的K更稳定。

我在项目中遇到过一件事：有个学生的数据只有12个样本，平行分析给出的K=1，但模型怎么调都不收敛。后来我把K手动设到3，问题直接解决了。原因就是小样本下平行分析不够敏感。

4.2 第二步：估计可替代变量

数量K定下来之后，接下来就是具体算这些变量了。SVA采用迭代加权SVD的方法，我把它叫做“剥洋葱”过程。

基本思路是：
先去除已知变量（比如实验条件）的影响，然后从残差里一步步提取出隐藏的结构。

我用伪代码来展示这个流程：

# 初始化
Y = 原始表达量矩阵 (n x m)
X = 已知协变量矩阵（包含我们关心的生物学变量）

# 第一步：回归掉已知协变量的影响
Y_resid = Y - X * (X^T X)^{-1} * X^T * Y

# 第二步：对残差矩阵做SVD
U, S, V = SVD(Y_resid)

# 第三步：提取前K个右奇异向量作为初始可替代变量
SV_initial = V[, 1:K]

# 第四步：迭代优化
repeat:
    # 将可替代变量加入模型，重新估计系数
    model = lm(Y ~ X + SV_initial)
    
    # 计算新的残差
    Y_new_resid = Y - fitted(model)[, !include_X_columns]
    
    # 对新残差做SVD
    U_new, S_new, V_new = SVD(Y_new_resid)
    
    # 更新可替代变量
    SV_new = V_new[, 1:K]
    
    # 检查收敛
    if ||SV_new - SV_initial|| < epsilon:
        break
    else:
        SV_initial = SV_new

# 输出最终的可替代变量
SV_final = SV_new

这个迭代过程的关键在于：
每一次迭代，都会把可替代变量和已知协变量一起回归掉，然后重新从残差中提取结构。这样就保证了可替代变量捕获的是真正独立的批次效应，而不是与已知变量混淆的混杂信号。

我记得第一次实现这个算法时，用了30次迭代没收敛。查了很久才发现是初始化SVD时用了原始数据而不是残差数据。这个小细节让我浪费了两天时间。

4.3 公式解读：为什么这样就能找到隐藏的批次效应？

如果你觉得上面的伪代码还不够直观，我们再看一下数学公式。

SVA假设的数据模型是：

Y_ij = μ_j + f_j(x_i) + ∑_{k=1}^{K} γ_jk * h_k(z_i) + ε_ij

其中：
- μ_j：基因j的基础表达水平
- f_j(x_i)：我们关心的生物学信号（比如疾病状态）
- γ_jk * h_k(z_i)：第k个可替代变量对基因j的影响
- ε_ij：随机噪声

SVA的目标就是：在不知道z_i的情况下，估计出h_k(z_i)（即可替代变量）和γ_jk（即基因负载）。

（1）当我们对残差矩阵做SVD时，实际上就是在对未知结构进行降维。
（2）SVD分解得到的右奇异向量，就是可替代变量的初始估计。
（3）迭代过程相当于在做“回归-残差-再回归”的交替最小二乘优化。

说白了，SVA就是把一个参数估计问题，拆解成了SVD+线性回归的组合。这个方法有一个明显的优势：完全不用假设可替代变量的分布形式，非常适合处理各种奇奇怪怪的批次效应。

4.4 本节小结

SVA的两步法设计得很精巧：
第一步用平行分析确定K，相当于先画个圈，知道要捕捉多少种隐藏效应。
第二步用迭代加权SVD，相当于把圈里的东西一个个捡出来。

我个人觉得，这种先粗后细的策略，比那些一次性估计K和变量的方法要稳健得多。至少在实际项目中，我用SVA处理过芯片数据、RNA-seq数据，甚至甲基化数据，K的估计都很稳定。