繰り返し測定した検査値をCox回帰へ入れる方法は、時変共変量で扱いました。start-stop形式に展開し、各区間に直近の測定値を置けば解析できます。
ただ、この方法では測定された値をそのまま共変量として使います。1回ごとの測定にばらつきがあっても、それを測定誤差としてモデルの中では分けません。また、次の測定までは前回値でその区間を代表させることになります。
血圧やCD4数のように測定ごとのばらつきがあり、時間とともに値も動く指標では、この違いが推定値に現れることがあります。
今回ランダムに生成したデータセットでは、直近の測定値を使ったCox回帰のハザード比が、測定誤差を大きくするにつれて1へ近づきました。結合モデルを使う理由は、この「観測値」と、その背後にある「患者ごとの軌跡」を分けて扱えるところにあります。
誤差が増えるとHRは0.56から0.82へ動いた
1,500例を最大5年間追跡し、3か月ごとに検査値を測定するデータをランダムに生成しました。イベントは443件、縦断測定は20,003回です。生存時間には観測値ではなく、その背後に設定した真の軌跡が関連するようにしてあります。
変更したのは測定誤差の標準偏差だけです。0から1.5まで動かし、その他の条件は同じにしました。
測定誤差がない場合のハザード比は0.56でした。誤差の標準偏差を1.0にすると0.74、1.5では0.82です。
0.56であれば、軌跡が1単位高い場合のハザードは44%低いという推定です。それが0.74になると26%低いという推定に変わり、対数ハザード比で見るとかなりの差であることが分かります。
説明変数の測定誤差によって関連が弱く推定される現象は、回帰の希薄化として知られています。ただし、ここは今回のシミュレーション結果と一般論を分けた方がよく、生存時間解析で時間依存性を持つ共変量を扱う場合、測定誤差があれば常に同じ方向・同じ大きさでハザード比が1へ寄る、とまではこの結果から言えません。
今回の設定では、測定間隔より測定誤差の影響の方が大きく出ました。
誤差の標準偏差を0.75に固定して測定間隔だけを変えると、ハザード比は1か月ごとで0.69、3か月で0.68、6か月で0.70、12か月でも0.70でした。測定回数を12倍にしても、ほとんど戻っていません。
この結果には、生成した軌跡が比較的緩やかだったことが効いています。区間内での真の変化より1回の測定に乗る誤差の方が大きい設定なら、測定間隔を短くしても誤差そのものは消えません。
逆に、短期間に大きく変化する検査値なら結果は違う可能性があります。したがって、「測定間隔を短くすれば希薄化を防げる」と一般化するのも、「測定間隔は関係ない」と一般化するのも避けた方が良いということが分かります。
結合モデルでは観測値ではなく軌跡を使う
検査値そのものではなく、測定誤差を除いた患者ごとの軌跡を生存時間と関連づけることが結合モデルの基本的な考え方です。
縦断データ側には、たとえば混合効果モデルを置きます。
$$y_i(t) = m_i(t) + \varepsilon_i(t), \qquad m_i(t) = \boldsymbol{\beta}^\top\mathbf{x}_i(t) + \mathbf{b}_i^\top\mathbf{z}_i(t) \tag{1}$$
$y_i(t)$は患者$i$の時点$t$で観測された値、$m_i(t)$はモデルから推定する潜在的な軌跡、$\varepsilon_i(t)$は測定ごとの誤差です。固定効果$\boldsymbol{\beta}$と患者ごとのランダム効果$\mathbf{b}_i$を使い、ばらついた測定値の背後にある軌跡を推定します。
生存時間側には、その$m_i(t)$を入れます。
$$h_i(t) = h_0(t)\exp\{\boldsymbol{\gamma}^\top\mathbf{w}_i + \alpha\, m_i(t)\} \tag{2}$$
$h_i(t)$は時点$t$のハザード、$h_0(t)$はベースラインハザード、$\mathbf{w}_i$はその他の共変量です。$\alpha$は、潜在軌跡$m_i(t)$とハザードの関連を表します。$\exp(\alpha)$を使えば、軌跡が1単位高い場合のハザード比として読めます。
式の形だけを見ると、時変Coxとそれほど違わないように見えますが、異なるのは右辺に入っている値です。
図2では、観測値が真の軌跡の上下に散らばっています。直近値をそのまま使えば、その散らばりも共変量として生存モデルへ入ります。
混合効果モデルでは、患者本人の複数回の測定と集団全体の情報を合わせて軌跡を推定します。そのため、1回だけ大きく外れた測定値の影響を、そのまま次回測定まで持ち越さずに済みます。
二段階法と同時推定は分けて考える
今回のランダムに生成したデータでは、推定軌跡をいったん求め、それをCox回帰へ渡す二段階法も計算しました。
ハザード比は、真の軌跡を使うと0.564、直近の観測値では0.679、混合効果モデルから推定した軌跡では0.582でした。
数字だけを見ると、二段階法でほとんど解決したように見えます。
ここで一つ問題があります。第1段階で得た$m_i(t)$は推定値ですが通常の二段階法では、それを第2段階のCox回帰で誤差のない既知の共変量のように扱います。第1段階の推定に伴う不確実性が、そのままでは第2段階の標準誤差に反映されません。
結合モデルは縦断サブモデルと生存サブモデルを同時に推定することで、この2つを一つの確率モデルとして扱います。WulfsohnとTsiatisが1997年に示したモデルも、この問題を扱った代表的な研究です。
私なら、解析を組む段階で二段階法を使うことはあります。軌跡を入れると点推定がどのくらい動くのかを短時間で確認できるからです。ただし、それをそのまま最終解析として採用するかは別に判断します。主要な推論対象が$\alpha$やその信頼区間なら、同時推定を優先します。
結び方も式(2)だけではありません。現在値$m_i(t)$を使うcurrent value associationのほか、変化速度を使う形、累積値を使う形、ランダム効果を生存サブモデルへ直接入れる形があります。
どの形を選ぶかで$\alpha$の意味は変わるため、値そのものとの関連を知りたいのか、低下速度との関連を知りたいのかを先に決めず、結果を見て最も都合のよい結び方を採用する進め方は推奨しません。
死亡で測定が終わるデータは欠測と少し違う
縦断データだけを見ていると、もう一つ問題があります。
状態の悪い患者ほど早くイベントを起こす研究では、その患者についてイベント後の縦断測定は得られません。死亡をイベントにしているなら、死亡後の血圧やCD4数は「本来は測れたが欠測した値」と同じものではありません。
そのため、通常の混合効果モデルだけで縦断データを解析し、死亡との関係を無視したまま平均軌跡を解釈するのは注意が必要です。
結合モデルでは、縦断過程とイベント時間を共有パラメータを介して一緒にモデル化します。死亡やイベントによって測定が終わる過程と縦断値の関連を、一定のモデル仮定のもとで扱えるところが利点です。
「死亡による脱落なら結合モデルを使えば自動的に補正される」と考えるのは適切ではなく、どのランダム効果を共有するか、縦断軌跡をどう書くか、生存側とどう結ぶかによって結果は変わります。
欠測データで扱った多重代入とも同じ問題ではありません。多重代入で通常扱う欠測と、終末イベント後には定義できない縦断アウトカムは分けて考えた方が整理しやすくなります。
動的予測は軌跡モデルの形に左右される
結合モデルには、ある時点までに得られた検査値の履歴を使い、その後のイベント確率を更新する使い方があります。Rizopoulosは2011年のBiometricsで、このdynamic predictionと予測精度の評価を扱っています。
今回の例では、2年時点まで生存している2人について、その時点までの縦断測定を使って予測を更新しました。
一方は2年時点の推定値が1.4で低下が速く、5年時点までの予測生存確率は0.16でした。もう一方は推定値5.4で低下が緩やかで、予測は0.94です。
登録時の検査値だけを固定して使うモデルとは違い、新しい測定が追加されれば、その情報を使って予測も更新できます。
ランドマーク解析も、ある時点まで生存した患者についてその後の予測を更新する方法です。結合モデルでは、そこまでに蓄積された縦断測定を軌跡としてモデル化できる点が違います。
ただし、予測が精密になるかどうかは結合モデルを使ったという事実だけでは決まりません。
たとえば実際の軌跡に強い非線形性があるのに、縦断サブモデルを直線で書けば、そのずれは動的予測にも入ります。スプラインや非線形な時間効果を使うべきかは、個々の観測値と推定軌跡を重ねて確認します。
予測目的なら、モデルを当てはめたデータで予測がよく見えるだけでは足りません。較正や時間依存的な予測精度を、内部検証または外部データで確認する必要があります。
結合モデルを使う理由を先に決める
繰り返し測定された変数があるからという理由だけで結合モデルを選ぶ必要はありません。
研究の中心が、縦断的に変化するバイオマーカーとイベント時間の関連にあるなら、結合モデルを使う理由は明確です。WulfsohnとTsiatisの測定誤差を含む縦断データと生存時間の研究や、その後の結合モデルの文献が扱ってきたのもこの問題です。
逆に、ある検査値を補助的な調整変数として入れたいだけで、測定誤差を明示的に扱うことが研究目的に直結しないなら、まず時変共変量を使ったモデルで十分かを考えます。ただし、「調整変数だから測定誤差を無視してよい」という意味ではありません。測定誤差によって主要な曝露の推定まで変わる可能性があるならば検討が必要です。
例えば、候補となる縦断マーカーが10個、20個とある状態で、片端から結合モデルを当てはめてp値で選ぶ方法も適切ではありません。どのマーカーを何の目的でモデル化するのかを先に絞る方が、最終的な$\alpha$の意味を説明できます。予測モデルとして多数の候補変数を扱う場合は、正則化Coxモデルで扱った変数数と過学習の問題も別に考える必要があります。
また、モデルが収束しない場合は、すぐにランダム効果を削る前に、時間とバイオマーカーの尺度を確認します。日単位の時間を年単位へ変える、値の桁が大きく違う変数を中心化・標準化する、といった変更だけで数値計算が安定することがあります。
加えて、測定回数も検討し、1回しか測定されていない患者から、その患者自身の時間的な傾きに関する情報は得られません。ただし、「ランダム傾きを使うには全患者が最低2回、安定推定には必ず3回必要」といった一律の基準があるわけではないため、1回測定の患者がどれくらいいるか、各時点に何人残っているか、ランダム効果の分散が推定できているかを合わせて判断します。
論文では縦断サブモデルと生存サブモデルの両方を書きます。縦断側は時間の関数形、固定効果、ランダム効果、測定誤差の仮定、生存側はベースラインハザードと共変量を記載します。そして、2つをcurrent value、slope、累積値などのどの形で結んだかを明記します。
観測値をそのまま使った時変Coxも併記すれば、今回の図1のように推定値がどの程度変わったかを読者が確認できます。ただし、値が結合モデルに近かったから時変Coxが正しく、離れたから間違いだった、という比較ではありません。推定しているモデルが違うことを前提に読ませる必要があります。
縦断サブモデルの当てはまりも図にします。何人かについて観測値と推定軌跡を重ねれば、直線で十分なのか、非線形な時間効果が必要なのかが分かります。モデルの形が最終的な推定値へ影響する点は、治癒率モデルで非治癒者の分布を変えたときとよく似ています。
実際の解析では、測定誤差の大きさ、測定回数、軌跡の形、生存側との結び方、そして主要な推論対象がどこにあるかを見て、時変Coxで止めるか結合モデルまで進むかを決めます。Dr.データサイエンスでは、こうした縦断データと生存時間を組み合わせた解析方針についてのご相談を承っています。
参考文献
- Wulfsohn MS, Tsiatis AA. A joint model for survival and longitudinal data measured with error. Biometrics. 1997;53(1):330-339.
- Rizopoulos D. Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data. Biometrics. 2011;67(3):819-829.
- Ibrahim JG, Chu H, Chen LM. Basic concepts and methods for joint models of longitudinal and survival data. Journal of Clinical Oncology. 2010;28(16):2796-2801.
- Rizopoulos D. Joint Models for Longitudinal and Time-to-Event Data: With Applications in R. Chapman and Hall/CRC; 2012.
- JMbayes2: Extended Joint Models for Longitudinal and Time-to-Event Data. CRAN.
