A scalable Bayesian functional factor model for high-dimensional longitudinal molecular data

这是一篇关于高维纵向分子数据可扩展贝叶斯功能因子模型的学术论文详细技术总结。该研究由苏黎世大学和剑桥大学医学研究委员会生物统计学小组的研究人员共同完成，旨在解决生物医学研究中大规模纵向分子谱分析的挑战。

以下是该论文的详细技术总结：

1. 研究背景与问题 (Problem)

• 背景：现代生物医学研究（如 COVID-19 纵向研究）产生了大量高维纵向分子数据（基因、蛋白质、代谢物等）。这些数据不仅维度极高（$p$ 可达数万），而且具有时间动态性。
• 核心挑战：
  1. 缺乏统一框架：现有的统计工具要么专注于纵向建模（如功能主成分分析 FPCA），要么专注于高维降维（如潜在因子模型），缺乏能同时处理变量间依赖和时间动态异质性的统一框架。
  2. 可解释性差：传统的多变量 FPCA 通常产生稠密的成分，难以解释具体的生物通路；而标准的因子模型往往忽略时间维度的平滑结构。
  3. 计算可扩展性：现有的基于高斯过程（Gaussian Process, GP）的因子模型（如 MEFISTO）在处理大规模样本（$N$）和高维变量（$p$）时，计算成本呈立方级增长，难以扩展。
  4. 数据稀疏性：临床数据通常采样稀疏且不规则，不同受试者的测量时间和次数不一致。

2. 方法论 (Methodology)

作者提出了一种可扩展的贝叶斯功能因子模型 (Scalable Bayesian Functional Factor Model)，其核心架构如下：

2.1 模型结构

模型结合了稀疏潜在因子建模与功能主成分分析 (FPCA)：

• 潜在因子表示：假设观测到的 $p$ 个功能变量 $y_i(t)$ 由少量 $Q$ 个潜在因子过程 $h_q(t)$ 线性组合而成：
  $$y_i(t) = \mu(t) + \sum_{q=1}^Q B_{jq} h_q(t) + \epsilon_i(t)$$
  其中 $B$ 是载荷矩阵，通过Spike-and-Slab 先验诱导稀疏性，确保每个潜在因子仅影响部分生物标志物（即特定的生物通路）。
• FPCA 分解：每个潜在因子过程 $h_q(t)$ 本身被建模为功能主成分分析的展开形式：
  $$h_q(t) = \sum_{\ell=1}^L \zeta_{i\ell}^{(q)} \psi_\ell^{(q)}(t)$$
  其中 $\psi_\ell^{(q)}(t)$ 是正交的特征函数（通过样条基表示），$\zeta_{i\ell}^{(q)}$ 是受试者特定的得分。这种分解将因子内的时间变异性分离为主要的正交模式。
• 贝叶斯设定：
  • 使用Spike-and-Slab 先验处理载荷矩阵 $B$ 的稀疏性，允许从数据中学习哪些变量属于哪些因子。
  • 使用半柯西先验 (Half-Cauchy) 处理方差参数，以增强鲁棒性。
  • 通过超参数 $\omega_q$ 控制每个因子的稀疏程度，并允许模型自动推断有效因子数量（通过过指定 $Q_{max}$ 并剔除不显著的因子）。

2.2 推断算法：退火变分推断 (Annealed Variational Inference)

由于后验分布具有多峰性（由 Spike-and-Slab 先验引起）且维度极高，传统的 MCMC 方法计算不可行。作者提出了一种退火变分推断 (Annealed VI) 算法：

• 变分推断 (VI)：使用平均场 (Mean-field) 假设将后验分布近似为可分解的分布，通过最大化证据下界 (ELBO) 进行优化。所有更新步骤均为闭式解 (Closed-form)，无需数值优化。
• 退火机制 (Simulated Annealing)：为了克服局部最优解问题，算法引入温度参数 $T \ge 1$。
  • 初始阶段设置较高的温度 ($T > 1$)，使目标函数表面变平，帮助算法探索多峰空间。
  • 随着迭代进行，温度逐渐降低至 $T=1$，最终收敛到原始变分目标。
  • 实验表明，退火机制显著提高了因子选择的稳定性和重构精度。
• 可扩展性：该算法避免了高斯过程协方差矩阵的求逆，计算复杂度随样本量 $N$ 和变量数 $p$ 近似线性增长，支持并行计算。

2.3 软件实现

该方法已实现为 R 包 bayesSYNC，支持并行执行，适用于处理数万变量和稀疏采样数据。

3. 主要贡献 (Key Contributions)

1. 统一的建模框架：首次将稀疏因子模型与功能主成分分析 (FPCA) 结合，能够同时捕捉跨变量的结构化依赖和因子内部的时间异质性。
2. 可解释性与稀疏性：通过 Spike-and-Slab 先验，模型能识别出特定的生物通路（即哪些标志物属于哪个因子），并量化个体在这些动态轴上的差异。
3. 可扩展的推断算法：开发了带有退火机制的变分推断算法，解决了高维、多峰后验分布的推断难题，使得在 $p=20,000$ 甚至更高维度下的分析成为可能。
4. 自动因子选择：利用过指定 (Overspecification) 策略和基于后验包含概率 (PPI) 的筛选机制，无需预先设定精确的因子数量即可自动识别有效因子。

4. 实验结果 (Results)

4.1 模拟研究

• 大规模数据表现：在 $N=100, p=20,000$ 的模拟数据中，模型准确恢复了潜在的时变结构、稀疏载荷模式以及有效因子数量。
• 退火机制的有效性：与标准变分推断 (Vanilla VB) 相比，退火变分推断显著提高了收敛质量，消除了无关因子的补偿效应，获得了更高的 ELBO 值和更准确的均值函数估计。
• 与 MEFISTO 的对比：
  • 准确性：在载荷恢复 (AUC) 和轨迹重构误差 (ISE) 上，bayesSYNC 优于基于高斯过程的 MEFISTO。
  • 可扩展性：MEFISTO 在处理 $N=40$ 时，部分运行因内存或时间限制失败（超过 36 小时），而 bayesSYNC 运行时间仅为几分钟，且随 $N$ 线性增长。
  • 可解释性：MEFISTO 默认假设所有受试者共享同一潜在轨迹（或仅通过高斯核耦合），无法像本模型那样通过 FPCA 分解出受试者特定的时间异质性模式。

4.2 应用：COVID-19 纵向分子数据

• 数据：剑桥医院 COVID-19 队列，包含 57 种炎症蛋白、细胞因子、代谢物等，覆盖症状出现后 7 周。
• 发现：
  • 因子 1 (炎症与色氨酸代谢)：识别出炎症标志物（CRP, IL-6 等）与色氨酸代谢通路（犬尿氨酸途径）的协同上调。这与长期新冠 (Long COVID) 中神经活性代谢物积累的理论一致。
  • 因子 2 (氨基酸耗竭)：识别出支链氨基酸（亮氨酸、异亮氨酸、缬氨酸）的耗竭模式，与重症患者的代谢特征相符。
  • 临床异质性：通过受试者层面的 FPCA 得分，模型量化了疾病严重程度和恢复速度的差异。得分高的患者表现出持续的炎症升高和恢复缓慢，而得分低的患者恢复较快。
  • 轨迹重构：模型成功重构了个体标志物的时间轨迹，并给出了 95% 可信区间，与临床观察高度一致。

5. 意义与结论 (Significance)

• 方法论突破：该研究提供了一种在超高维、稀疏纵向数据中进行贝叶斯推断的可行方案，解决了传统高斯过程模型在计算上的瓶颈。
• 生物医学价值：模型不仅能降维，还能提供具有生物学意义的解释（如通路特异性），帮助研究人员理解疾病进展中的动态机制和个体差异。
• 临床转化潜力：生成的低维受试者得分（Subject-level scores）可作为生物标志物，用于预测疾病严重程度、监测治疗反应或识别长期后遗症风险。
• 未来方向：作者指出，未来可将协变量直接纳入因子得分分布以关联临床结果，并进一步探索高维功能设置下的后验理论性质。

总结：这篇论文通过结合稀疏因子模型、FPCA 和退火变分推断，成功构建了一个可扩展、可解释且计算高效的框架，为分析现代生物医学中日益复杂的高维纵向分子数据提供了强有力的工具。