後ろ向きのデータベース研究で、ある治療を受けた患者と受けなかった患者の生存を比べるケースがあります。治療群のほうが明らかに長く生きている。ハザード比は0.3を切り、信頼区間も1をまたがない。
ここで確かめておきたいことがあります。その治療を受けるまでに、どのくらいの時間がかかっていたか。移植なら待機期間、二次治療なら一次治療の期間、薬剤なら処方が始まるまでの期間です。その間に亡くなった患者は、定義上「治療を受けなかった群」に入ります。
つまり治療群には、治療を受けるまで生き延びたという条件が最初から埋め込まれています。この期間を不死時間と呼びます。
待機している間は死ねない
この構造がどれくらいの誤差を生むかは、実データで確かめられます。survivalパッケージのjasaは1967年から1974年のスタンフォード心移植プログラムの記録で、この問題が最初に指摘されたときのデータそのものです。103例のうち69例が移植を受け、待機日数の中央値は25日、最長は309日でした。
library(survival)
candidate <- jasa
candidate$transplanted <- as.numeric(!is.na(candidate$tx.date))
candidate$followup <- as.numeric(candidate$fu.date - candidate$accept.dt)
coxph(Surv(followup, fustat) ~ transplanted, data = candidate)#> exp(coef) lower .95 upper .95
#> transplanted 0.266 0.165 0.429ハザード比0.27。移植を受けると死亡ハザードが7割下がる、という結果です。
ところが移植群69例の待機日数を合計すると2600日あり、この時間は全て移植群の観察時間として数えられています。しかも移植を受けなかった34例の追跡日数の中央値は20日で、多くが待機中に亡くなっています。移植群と非移植群の違いには、移植の効果と「25日生き延びられた」という選択の両方が混ざっています。
共変量を時間の関数にする
解決は、移植の有無を追跡開始時点の属性ではなく、時点ごとに変わる状態として扱うことです。Coxモデルの共変量を時間の関数に置き換えます。
$$h(t \mid \mathbf{Z}(t)) = h_0(t)\exp(\boldsymbol{\beta}^\top \mathbf{Z}(t)) \tag{1}$$
式(1)では、同じ患者が時点によって違う共変量値を持ちます。移植前の期間は$Z(t)=0$、移植後は$Z(t)=1$です。部分尤度の各因子は、その時刻のリスク集合にいる人の「その時点での」共変量値で計算されます。
$$L_j(\boldsymbol{\beta}) = \frac{\exp(\boldsymbol{\beta}^\top \mathbf{Z}_{(j)}(t_j))}{\sum_{i \in R(t_j)} \exp(\boldsymbol{\beta}^\top \mathbf{Z}_i(t_j))} \tag{2}$$
式(2)が意味するのは、待機中の患者は「まだ移植を受けていない人」として分母に入るということです。移植前に亡くなれば、その死亡は移植前の状態に帰属します。不死時間が移植群に付け替えられることがなくなります。
データの持ち方も変わります。1人1行ではなく、状態が変わるたびに行を分けます。
この形をstart-stop形式と呼びます。Surv関数に開始時点と終了時点の両方を渡すと、Rはこの構造を読み取ります。左切断のときに使った3引数の書き方と同じです。
head(heart[, c("id", "start", "stop", "event", "transplant")], 4)#> id start stop event transplant
#> 1 1 0 50 1 0
#> 2 2 0 6 1 0
#> 3 3 0 1 0 0
#> 4 3 1 16 1 1coxph(Surv(start, stop, event) ~ transplant, data = heart)#> exp(coef) lower .95 upper .95 Pr(>|z|)
#> transplant1 1.136 0.629 2.049 0.673
0.27が1.14になりました。信頼区間も0.63から2.05で1を含みます。同じデータ、同じイベント数で、結論が逆になります。差はデータの持ち方だけです。
この解析は1970年代前半に議論され、Crowley and Hu（J Am Stat Assoc 1977）が時変共変量による扱いを示しました。50年前に決着した話ですが、同じ誤りは今も繰り返されています。Suissa（Am J Epidemiol 2008）は、薬剤疫学の観察研究で不死時間バイアスがどのように混入するかを整理しています。
検査値を入れるときは話が別
移植の有無のように、患者の状態とは別のところで決まる共変量は外因性と呼ばれます。治療の切り替え、季節、環境曝露。これらは時変共変量として素直に扱えます。
一方、CD4数、腫瘍マーカー、eGFR、症状スコアのように、疾患そのものが生み出す値は内因性です。ここに同じ手続きを当てはめると、別の問題が出ます。
時点$t$で測定された検査値は、その患者が$t$まで生存したという条件のもとで観測されています。悪化した値を持つ患者は、そもそも測定される前に亡くなっていることがあります。ハザード比を推定しても、それは「検査値が悪いと死亡が近い」という関連を表すだけで、検査値を改善させれば死亡が減るという意味にはなりません。
さらに厄介なのは、治療効果を見たいときに内因性の検査値を調整に入れる場合です。治療が検査値を通じて生存を改善しているなら、その検査値は治療と転帰の中間にあります。中間変数で調整すると治療効果の一部が消えるので、過調整になります。
そのため、内因性の時変共変量をCoxモデルに入れて得たハザード比を、そのまま因果的に読むことは勧めません。予後の記述として使うか、あるいは縦断データと生存データの結合モデルや逆確率重み付けのような、この構造を明示的に扱う手法に進むことになります。
もうひとつ、測定のタイミングにも注意が要ります。状態が悪くなった患者ほど検査の頻度が上がるなら、測定の有無自体が予後の情報を持ちます。区間内で共変量を一定と見なす近似も、測定間隔が長いほど粗くなります。
ランドマーク解析という別の道
時変共変量を使わずに不死時間を避ける方法もあります。基準時点$t_L$を決め、そこまで生存した患者だけを対象にして、$t_L$時点での状態で群分けする。ランドマーク解析です。
移植の例なら、登録から60日を基準時点にして、60日時点で生存していた患者を、そこまでに移植を受けたかどうかで分けます。60日より前に亡くなった患者は両群から除かれるので、不死時間は生じません。
単純で、読者にも説明しやすい方法です。ただし$t_L$より前に起きたイベントの情報を捨てますし、$t_L$をどこに置くかで結果が変わります。私なら、$t_L$を臨床的な根拠（治療の標準的な導入時期など）から決めたうえで、複数の$t_L$で感度分析を出します。結果を見てから$t_L$を選ぶと、変数選択と同じ問題が起きます。
状態が変わった日付が分からないとき
時変共変量の解析には、状態が変わった日付が要ります。この日付がデータにないことが、実務では頻繁に起きます。
後ろ向きの抽出では、現在の状態しか出てこないことがあります。「この患者は二次治療を受けた」という情報はあるが、いつ移行したかがない。この状態では時変共変量として扱えません。
代替として使われるのがランドマーク解析です。基準時点での状態だけを使うので、変化日が不要になります。ただし基準時点より前の変化と後の変化を区別できないので、情報は落ちます。
もうひとつ、日付が一部の患者でしか分からない場合もあります。分かる患者だけで解析すると、記録が残っている患者に偏ります。記録の有無が診療の丁寧さと関係していれば、選択バイアスになります。この場合は、欠測の割合と、欠測した患者の特徴を報告することになります。
区間を細かく切りすぎない
start-stop形式では、区間の切り方を解析者が決めます。共変量が変わるたびに切るのが基本ですが、連続的に変わる変数では細かく切りたくなります。
細かくすると行数が増えます。1人あたり100区間なら、1000例で10万行です。計算時間とメモリを消費しますが、推定値はほとんど変わりません。部分尤度が使うのはイベント時刻でのリスク集合なので、イベントが起きない区間をいくら細かく切っても情報は増えないためです。
実務的には、共変量が実際に変化した時点でのみ切れば足ります。検査値なら測定日、治療なら開始日と終了日です。等間隔で機械的に切る必要はありません。
データを受け取ったときに確認すること
曝露や治療の状態が追跡期間中に変わるなら、それがいつ変わったかの日付が必要です。この日付がないデータでは、時変共変量の解析はできません。後ろ向きの抽出を依頼する段階で、状態の変化日を含めてもらうことになります。
1人1行のデータを受け取ったとき、曝露の定義に「追跡期間中に一度でも」という言葉が入っていたら、不死時間バイアスを疑う場面です。「入院中にこの薬剤を投与された患者」「経過中に二次治療へ移行した患者」といった分類がこれにあたります。
報告では、start-stop形式に展開したこと、状態の変化をどの日付で定義したか、そして各群の人数ではなく人年を書きます。時変共変量では患者が両群に登場するので、人数だけでは追跡の実態が伝わりません。
比例ハザード性の確認は時変共変量でも必要で、扱いは比例ハザード仮定の検証と同じです。モデルの基本形はCox比例ハザードモデルにあります。
後ろ向きのデータで曝露の定義をどう置くか、ランドマーク解析と時変共変量のどちらを主解析にするか。このあたりはデータの取得可能性と研究目的の両方を見ないと決まりません。Dr.データサイエンスでは、こうした解析設計のご相談を承っています。
参考文献
Crowley J, Hu M. Covariance analysis of heart transplant survival data. Journal of the American Statistical Association. 1977;72(357):27-36.
Gail MH. Does cardiac transplantation prolong life? A reassessment. Annals of Internal Medicine. 1972;76(5):815-817.
Suissa S. Immortal time bias in pharmacoepidemiology. American Journal of Epidemiology. 2008;167(4):492-499.
Dafni U. Landmark analysis at the 25-year landmark point. Circulation: Cardiovascular Quality and Outcomes. 2011;4(3):363-371.
Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer; 2000.
