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

1. 問題定義と背景

現代の生物医学研究では、個体ごとに遺伝子、タンパク質、代謝物などの分子マーカーを時間経過とともに繰り返し測定する「大規模な縦断分子プロファイリング」が一般的になっています。しかし、これらのデータを分析するには以下の課題が存在します。

• 高次元性と複雑な依存関係: 変数（分子マーカー）の数が非常に多く（数万件規模）、関連する生物学的経路内のマーカー間には強い相関があります。
• 時間的ダイナミクスと個人差: 生物学的プロセスは時間とともに変化し、患者間ではその発現パターン（回復速度や重症度）に大きな個人差があります。
• 既存手法の限界:
  • 潜在因子モデル: 高次元データの次元削減に有効ですが、通常は観測値を多変量として扱い、明示的な「関数（時間的連続性）」としてモデル化しないため、時間的ダイナミクスを捉えきれません。
  • 関数主成分分析（FPCA）: 時間的変動を捉えますが、通常は変数ごとに独立して分析するか、全変数に負荷を持つ「密（dense）」な成分を生成するため、生物学的な解釈（どのマーカーがどの経路に関与しているか）が困難です。
  • MEFISTO（ガウス過程ベース）: 柔軟な時間モデルを提供しますが、計算コストが高く、大規模なサンプル数や変数数に対してスケーラビリティに課題があります。また、変数ごとのスパースな構造を明示的にモデル化しないため、解釈性が低い場合があります。

目的: 変数間の共依存関係と、プロセス内の時間的異質性を同時にモデル化し、かつ数千〜数万の変数に対しても計算的に実行可能で、生物学的に解釈可能な枠組みの構築。

2. 提案手法：ベイズ関数因子モデル

提案されたモデルは、スパースな潜在因子モデルと**関数主成分分析（FPCA）**を階層的に結合したものです。

モデルの構造

1. 関数因子表現:
   観測された $p$ 次元の関数データ $y_i(t)$ を、低次元の潜在因子プロセス $h_q(t)$ の線形結合として表現します。
   $$y_i(t) = \mu(t) + \sum_{q=1}^Q b_{jq} h_q(t) + \epsilon_i(t)$$
   ここで、$b_{jq}$ は因子 $q$ と変数 $j$ の間の**負荷（loading）**です。

2. スパース性の導入（Spike-and-Slab Prior）:
   負荷行列 $B$ に対してスパース性を誘導するため、スパイク・アンド・スラブ事前分布を採用します。これにより、各因子が関連する変数のサブセットのみを選択的に負荷し、生物学的経路ごとの明確な構造を抽出します。
   $$b_{jq} | \gamma_{jq} \sim \gamma_{jq} \mathcal{N}(0, 1) + (1-\gamma_{jq})\delta_0$$

3. FPCA による時間的変動の分解:
   各潜在因子プロセス $h_q(t)$ 自体も時間的に変動するため、これを FPCA を用いて直交モードに分解します。
   $$h_q(t) = \sum_{\ell=1}^L \zeta_{i\ell}^{(q)} \psi_\ell^{(q)}(t)$$
   ここで、$\psi_\ell^{(q)}(t)$ は固有関数（時間的パターン）、$\zeta_{i\ell}^{(q)}$ は個体 $i$ のスコア（そのパターンへの個人差の度合い）です。これにより、因子ごとの主要な時間的パターンと、個体ごとの変動を分離して記述できます。

4. スケーラブルな推論アルゴリズム（Annealed Variational Inference）:
   高次元かつスパースな事後分布の推論は計算的に困難です。そこで、**変分推論（Variational Inference, VI）を採用し、さらにシミュレーテッド・アニーリング（Annealing）**を組み合わせたアルゴリズムを開発しました。

   • 変分推論: 事後分布を近似する分布を最適化することで、MCMC に比べて高速な推論を可能にします。
   • アニーリング: 負荷行列のスパース構造により生じる多峰性の事後分布において、局所最適解に陥るのを防ぐため、温度パラメータを徐々に下げながら最適化を行います。これにより、より安定した因子選択と精度向上を実現しています。

3. 主要な貢献

1. 統合的な枠組みの提案: 高次元縦断データに対して、変数間のスパースな依存関係と、時間的ダイナミクス（FPCA）を統一的にモデル化する初めての手法の一つです。
2. スケーラビリティの実現: 変数数 $p$ が 20,000 に達するシミュレーションデータでも正確に推論可能なアルゴリズムを開発しました。従来のガウス過程ベースの手法（MEFISTO など）は、サンプル数が増えると計算コストが爆発的に増大しますが、本手法は線形に近いスケーリングを示します。
3. 解釈性の向上:
   • スパースな負荷: どの分子マーカーがどの生物学的経路（因子）に関与しているかを明確に特定できます。
   • FPCA スコア: 個体の時間的変動パターンを低次元のスコアで要約し、臨床的な重症度や回復動態との関連を定量化できます。
4. ソフトウェアの公開: 手法を R パッケージ bayesSYNC として公開し、実用性を高めています。

4. 結果

シミュレーション研究

• 大規模データでの性能: $N=100$ 個体、$p=20,000$ 変数のデータにおいて、真の因子数、スパースな負荷パターン、時間的構造（固有関数）を高精度に復元しました。
• アニーリングの重要性: 通常の（アニーリングなしの）変分推論と比較し、アニーリングを導入することで、不要な因子の選択を防ぎ、推定精度（ISE）と ELBO（証拠下限）が向上することを実証しました。
• MEFISTO との比較:
  • 精度: 負荷の復元精度（AUC）や軌道の再構成誤差（ISE）において、MEFISTO よりも優れた性能を示しました。特に、スパース性を明示的にモデル化している点が有利に働いています。
  • 計算コスト: サンプル数 $N$ が増加するにつれ、MEFISTO（個体ごとの軌道を推定する設定）は計算時間が急増し、メモリ不足や時間制限に直面しましたが、本手法は $N=40$ でも数分以内で完了しました。

COVID-19 縦断データへの適用

Cambridge 大学病院で行われた COVID-19 患者の分子プロファイリングデータ（57 種類のバイオマーカー、7 週間の追跡）に適用しました。

• 因子の発見: 2 つの主要な潜在因子が抽出されました。
  • 因子 1（炎症・キヌレニン経路）: CRP、サイトカイン、補体成分、およびキヌレニン経路代謝物（キヌレニン、キノリン酸など）と正の相関、トリプトファンやセロトニンと負の相関を示しました。これは、炎症性サイトカインによる IDO 酵素の活性化とトリプトファン分解経路の活性化という既知の生物学的メカニズムと一致します。
  • 因子 2（アミノ酸代謝）: 分岐鎖アミノ酸（ロイシン、イソロイシン、バリン）など、アミノ酸代謝経路に関連するマーカーが支配的でした。
• 臨床的解釈:
  • 因子 1 の固有関数の第 1 成分は「重症度」の指標となり、重症患者ほど炎症マーカーの持続的な上昇を示しました。
  • 第 2 成分は「回復動態」を捉え、スコアが大きい患者ほど早期にベースラインへ回帰しました。
  • これらのスコアは、患者の臨床分類（無症状、軽症、重症、人工呼吸器使用など）と有意に相関しており、個体間の回復プロセスの異質性を定量的に記述できました。

5. 意義と結論

この研究は、大規模な縦断分子データ解析において、**「計算効率性」「統計的精度」「生物学的解釈性」**のバランスを達成する新しいパラダイムを提供しました。

• 方法論的意義: ガウス過程ベースの柔軟性と、因子モデルのスパース性・解釈性を組み合わせ、変分推論とアニーリングによって大規模データへの適用を可能にしました。
• 応用可能性: 自己免疫疾患、がんの進行、治療反応、発達プロセスなど、時間的組織化が重要なあらゆる生物医学的縦断研究に応用可能です。
• 臨床的インパクト: 個体ごとの回復動態を低次元スコアで要約することで、個別化医療や予後予測への応用が期待されます。

結論として、この手法は複雑な生物学的メカニズムを解明し、患者の異質性を理解するための強力なツールとなり得ます。