論文の技術的サマリー：縦断データと時間至イベントデータの共同モデリングにおける SNP 全体効果の評価のためのサンプルサイズと検出力決定

1. 研究の背景と問題提起

臨床研究、特に糖尿病などの慢性疾患の遺伝子研究において、ヘモグロビン A1c（HbA1c）のような縦断的バイオマーカーは、時間至イベント（生存）アウトカムと強く関連しているため頻繁に収集されます。これらのデータを統合して分析するための「縦断データと生存データの共同モデリング（Joint Modeling）」手法は既に確立されています。

しかし、遺伝子研究における統計的デザイン、特にサンプルサイズと検出力の計算に関する検討は、既存の生存分析や治療効果の評価に比べて十分ではありません。遺伝子変異（SNP）は、イベントリスクに直接影響を与えるだけでなく、縦断的バイオマーカーを介した間接的な影響も及ぼす可能性があります。従来の方法では、SNP のようなカテゴリカル変数（0, 1, 2 コピー）や、両側検定が必要な遺伝的効果の特性を適切に反映したサンプルサイズ計算式が不足していました。

2. 提案手法と方法論

2.1 統計的枠組み

本研究では、縦断プロセスと生存プロセスをリンクさせる共同モデリング枠組みを採用しました。

• 生存モデル: 比例ハザードモデルを用い、ハザード関数に真の縦断軌道（latent trajectory）を含めます。
  $$\lambda_i(t) = \lambda_0(t) \exp \{ \gamma_g \text{SNP}_i + \alpha \eta_i(t) \}$$
  ここで、$\gamma_g$ は SNP の直接効果、$\alpha$ は縦断軌道とイベントリスクの関連を表します。
• 縦断モデル: 個体ごとの多項式関数で真の軌道 $\eta_i(t)$ をモデル化し、SNP の影響（$\beta_g$）を含めます。
  $$\eta_i(t) = (\beta_0 + b_{i0}) + (\beta_1 + b_{i1})t + \dots + \beta_g \text{SNP}_i + \epsilon_i(t)$$

2.2 全体効果の定義

SNP の生存アウトカムへの全体効果は、直接効果と間接効果（縦断プロセスを介した効果）の和として定義されます：
$$\text{全体効果} = \gamma_g + \alpha \beta_g$$
本研究の主要な目的は、この全体効果（$\gamma_g + \alpha \beta_g = 0$）を検証するための検定を行う際のサンプルサイズを決定することです。

2.3 サンプルサイズ計算式の導出

Schoenfeld (1983) や Chen et al. (2011) のアプローチを拡張し、SNP の特性（3 値のカテゴリカル変数、両側検定）を考慮した閉形式（closed-form）のサンプルサイズ計算式を導出しました。

必要なイベント数 $D$ は以下の式で与えられます：
$$D = \frac{(Z_{1-\tilde{\beta}} + Z_{1-\tilde{\alpha}/2})^2}{2pq(\gamma_g + \alpha\beta_g)^2}$$

• $p, q$: 対立遺伝子頻度（$q=1-p$）。SNP の分散 $2pq$ が式に含まれる点が、二値変数の治療効果の式との主な違いです。
• $\tilde{\alpha}$: 有意水準（両側検定のため $\tilde{\alpha}/2$ を使用）。
• $\tilde{\beta}$: 検出力（$1-\tilde{\beta}$）。
• $\gamma_g + \alpha\beta_g$: SNP の全体効果サイズ。

この式は、対立遺伝子頻度が低い（$p$ が小さい）場合や、効果サイズが小さい場合、より大きなサンプルサイズが必要になることを示しています。

2.4 実用的な推定アプローチ（2 ステップ法）

共同モデリングの直接推定は計算コストが高く、特にゲノムワイド解析では現実的ではありません。そのため、本研究では計算負荷を軽減する2 ステップ法を提案し、検出力評価に利用しました。

1. 縦断モデルに SNP を含めて $\hat{\beta}_g$ を推定し、観測値から SNP 関連の変動を除去した調整済み応答 $y'_i(t)$ を作成する。
2. 調整済み軌道 $\hat{\eta}'_i(t)$ を用いて生存モデルをフィットし、全体効果 $\theta_g$（$\gamma_g + \alpha\beta_g$ の推定値）を検定する。

3. シミュレーション研究と結果

3.1 計算式の精度と頑健性

1000 回のシミュレーションを行い、提案された計算式が有限サンプルにおいて正確かつ頑健であることを確認しました。

• パラメータ変化への耐性: ランダム効果の分散構造（$\delta^2_0, \delta^2_1$）や、$\alpha, \beta_g, \gamma_g$ の値が異なっても、計算された検出力とシミュレーションによる実測検出力はよく一致しました。
• 測定頻度: 十分な測定回数がある限り、測定スケジュール（頻度）の違いには比較的 insensitive でした。ただし、測定回数が極めて少ない場合（区間打ち切りに近い状況）、検出力の過小評価が生じる傾向がありました。
• モデルの誤指定: 真の縦断軌道が二次関数である場合、線形モデルを誤って指定すると、実測検出力が計算値から乖離し、低下することが示されました。これは、モデルの誤指定が検出力推定の精度に影響を与えることを示唆しています。

3.2 バイアスの評価

縦断プロセスを生存モデルから除外した場合（Cox モデルのみ）、$\gamma_g$ の推定にバイアスが生じることが確認されました。特に、SNP が縦断プロセスに影響を与え（$\beta_g \neq 0$）、かつ縦断プロセスが生存に影響を与える（$\alpha \neq 0$）場合、そのバイアスは顕著になります。提案する 2 ステップ法や共同尤度法は、このバイアスを適切に制御できることを示しました。

3.3 パラメータと検出力の関係

• 追跡期間: 追跡期間が長いほど、イベント数が増加し、検出力は向上しました。
• 効果サイズと対立遺伝子頻度: 全体効果サイズ（$\gamma_g + \alpha\beta_g$）が小さい、または対立遺伝子頻度が低い（稀な変異）場合、同等の検出力を得るためには大幅にサンプルサイズを増やす必要があります。

4. 実データへの適用：DCCT 研究

提案手法の実用性を検証するため、糖尿病管理と合併症試験（DCCT）のデータを用いた回顧的検出力分析を行いました。

• データ概要: 1 型糖尿病患者 1,441 名（集約治療群 527 名、標準治療群 526 名）のデータ。HbA1c の縦断測定と網膜症の発症時間を分析対象としました。
• 結果: 従来の GWAS 有意水準（例：$5 \times 10^{-8}$）や標準的な有意水準（0.05）において、DCCT のサンプルサイズでは、中等度の遺伝的効果を検出するための十分な検出力が得られていないことが示されました。
  • 標準治療群では、対立遺伝子頻度 $p > 0.20$ で 80% の検出力が得られますが、集約治療群（イベント数が少ない）では $p \approx 0.40$ が必要でした。
  • この結果は、既存の DCCT 解析で検出された遺伝的関連が、統計的検出力の限界により見逃されている可能性や、より大規模なコホートが必要であることを示しています。

5. 結論と意義

主要な貢献

1. 閉形式の計算式の提供: 縦断データと生存データを共同モデル化する枠組みにおいて、SNP の全体効果（直接＋間接）を検定するための、実用的なサンプルサイズ・検出力計算式を初めて導出しました。
2. SNP 特有の課題への対応: カテゴリカル変数（3 値）と両側検定を考慮し、対立遺伝子頻度を明示的に組み込んだ式を提供しました。
3. ツールの開発: 計算式を実装した対話型 Shiny アプリ（PowerSNP）を開発し、研究者が容易に利用できるようにしました。

意義と限界

• 意義: 遺伝子研究の計画段階において、縦断的バイオマーカーを介した間接効果を含む SNP の検出に必要なサンプルサイズを事前に正確に見積もることを可能にし、研究デザインの効率化と統計的検出力の向上に寄与します。
• 限界と注意点:
  • 式は「非情報的中止（non-informative censoring）」を仮定しており、情報的中止がある場合はバイアスのリスクがあります。
  • 縦断モデルの誤指定（例：真の非線形性を線形モデルで近似するなど）は、パラメータ推定と検出力計算の精度に影響を与えます。
  • 間接効果と直接効果が逆方向に働く場合、解釈や検出力推定が複雑になる可能性があります。

本研究は、複雑な遺伝子 - 環境相互作用を含む縦断・生存データ分析において、統計的デザインの厳密性を高めるための重要な基盤を提供しています。