📚 课程目录第3章 MOFA模型参数

MOFA 模型参数：因子数量、协方差结构与先验分布

说实话，多因子分析里最让人头疼的，就是那几个核心参数怎么设。许多人一上来就调，结果模型要么不收敛，要么解释不通。我记得有一次，一个项目跑了两周才发现因子数设太多了，回头重来……哎，真是血泪教训。

这一章，我就把自己踩过的坑、总结出来的经验，都摊开来讲。核心就三个东西：因子数量（K）、协方差结构和先验分布。这三者决定了你的模型是“靠谱的生物学发现”还是“一堆噪音”。

因子数量（K）—— 选多了还是选少了？

因子数量说白了，就是你推测数据背后藏着多少个隐藏的“基架”模式。它直接影响模型的复杂度和解释力。

K 值的确定策略

我个人习惯，先用“总方差解释比例”来做快速筛查。设定一个阈值，比如 80%~90%，看需要多少因子才能达到这个比例。

# 快速评估因子数量（伪代码示意）
threshold = 0.85
cumulative_variance = 0.0
k = 0
while cumulative_variance < threshold:
    k += 1
    cumulative_variance += explained_variance[k-1]
print(f"建议因子数: {k}")

不过这只能给个大致范围。更严谨的做法，是用坦普尔-韦斯特统计量（TWW）。它衡量的是“增加一个因子后，模型拟合改善了多少”，还考虑了你惩罚了模型的复杂度。嗯，这个指标在 MOFA2 的 Python 包里直接就提供了。

实用小技巧：我一般会从K=5开始扫，一直扫到K=30或者样本数/5（取小的那个）。然后画出 TWW 曲线，选那个“拐点”位置——也就是再增加因子数，收益明显下降的地方。

交叉验证法

如果数据量够，我会做 k-fold 交叉验证。把数据分成几份，留一份当验证集，剩下的训练模型。然后看预测误差。

# 交叉验证示意
from mofa import MOFA2

best_k = None
lowest_error = float('inf')
for k in candidate_ks:
    model = MOFA2(n_factors=k)
    model.fit(train_data)
    error = model.evaluate(val_data)
    if error < lowest_error:
        lowest_error = error
        best_k = k

避坑指南：我曾经在一次蛋白质组+转录组的多组学整合项目里，把K设到了25。结果跑出来十几个“因子”解释的全是批次效应。后来用交叉验证大幅削减到8个因子，才真正抓到生物学信号。所以不要贪多，K值宁少勿滥。

协方差结构 —— 给因子之间的“关系”上规矩

协方差结构定义的是因子的载荷（loading）之间，允不允许有关系，有什么关系。

四种常见结构

结构类型	描述	自由度	适用场景
零矩阵	各特征的载荷之间互相独立	0	特征间无共享变异
单位矩阵	所有因子载荷的方差相等且独立	0	简单基准模型
对角矩阵	不同特征的载荷方差不同，但彼此独立	特征数	允许特征有不同噪音水平
满秩矩阵	载荷之间可以相关	特征数×(特征数+1)/2	特征间存在复杂的相关结构

你想想看，“满秩”听起来很厉害对吧？但实际上它参数太多。没个几千样本，根本估计不出来。我一般会这么选：

数据量小（样本 < 50）：直接上零矩阵或对角矩阵，避免过拟合。
中等数据（样本 50-200）：用对角矩阵，如果是单细胞数据甚至可以手动设定对角加一个小的约束项。
大样本（样本 > 200）：可以考虑满秩，但还是要配合交叉验证看泛化能力。

一句话总结：协方差结构越复杂，需要的样本量越大。别为了一点点拟合优度，把模型弄塌了。

先验分布 —— 告诉贝叶斯“你该信什么”

MOFA 是贝叶斯模型，所以每个参数都要有一个先验。先验说白了，是我们在看数据之前，对这个世界已有的看法。但它最容易被忽略，却最能出问题。

核心参数先验

MOFA 里的先验主要落在“方差解释量”和“载荷的稀疏程度”上。关键参数有两个：

方差阈值（variance_threshold）：控制特征被至少一个因子解释的最小方差比例。默认是 0.1%，但我经常调到 1% 甚至 5%——特别是在做了大量的差异分析之后。
自由度（dof）：控制载荷的先验分布形状。dof 越小，先验越厚尾，越容易产生很大的估计值。

我常用的先验设定

说实话，我自己摸索出来的套路是这样的：

# MOFA2 Python 包中的先验设置示例
from mofax import MOFA

model = MOFA()
model.set_data(data)
model.set_model_options(
    n_factors=10,
    spikeslab_weights=True,  # 使用 spike-and-slab 先验来做变量选择
    ard_weights=True,        # 自动相关性决定，做因子稀疏化
    ard_variance=True,       # 对方差做自动相关性决定
    variance_threshold=0.01  # 每个因子至少解释 1% 的方差
)
model.build_model()

这里有两层先验：

Spike-and-slab 先验：让每个因子对大多数特征的载荷趋向于0，只保留少数特征的强关联。说白了，就是让因子解读起来更清晰。
ARD 先验：自动相关性决定。它会主动把不重要的因子“推”向0，从而帮你自动收缩因子数量。

个人经验：我通常会把 ard_weights=True 打开，因为它能帮你自动做变量选择。但注意了，如果数据量很小（比如几十个样本），ARD 可能会过于激进地把所有因子都压缩掉。这时候我会把 ard_variance 先关掉，或者给一个更大的先验规模。

参数调优的总思路

讲到这里，你可能发现这三者其实是联动的。因子数量多，协方差结构就不能太复杂；先验设得严，因子数就会自动收敛。下面这张图，是我自己整理的选参数流程：

这个流程里最关键的，是先验设定那一步。别小看它，我曾经见过有人把 ard_weights 设错了，结果模型识别出 30 个“显著”因子，但验证集上一塌糊涂。后来发现是因为 ARD 的先验尺度设得太大，所有因子都被判为“有贡献”。

一个真实翻车案例：在一次 RNA-seq + miRNA-seq 的整合分析中，我最初用了默认的零矩阵协方差结构。结果两个组学的因子载荷完全不相关，生物功能根本对不上。后来换成对角矩阵，并给 miRNA 的载荷加了更弱的先验（更大的方差），才得到一组在两条通路共有的因子。所以你想想看，协方差和先验是联动的，不能孤立调。

写在参数选择之后

参数这个东西，没有绝对正确的。我个人习惯是先跑一个简单的基准模型（K=10，零矩阵，默认先验），看收敛性和基础方差解释比例。然后一步步放松约束，每次只动一个参数，用交叉验证打分。这样虽然慢，但能保证每一步改动都有依据。

最终，别忘了再回到你的生物学假设里去验证。MOFA 给的因子，最好能和已知的通路、细胞类型或实验条件对得上。如果对不上，八成是参数没选对——而不是生物学错了。