レセプトや電子カルテを使った研究で、ある薬を使った患者と使わなかった患者の生存を比べることがあります。追跡期間のどこかで処方があれば治療群、なければ対照群という分け方は、データの持ち方として自然に見えます。ところが、この分け方で5000例の架空データを解析すると、ハザード比は0.42になりました。データを作った時に置いた真の効果は0.70です。効果が実際の2倍近く大きく見えていることになります。
ズレの原因は効果の推定方法ではなく、追跡の起点をどこに置いたかにあります。
時変共変量で扱った不死時間バイアスが、データベース研究では解析の技法ではなく設計そのものの問題として現れます。
治療群に入るには生き延びる必要がある
登録から8か月後に治療を始めた患者を考えます。この患者が治療群に入るためには、8か月間は死なずにいなければなりません。一方、3か月で亡くなった患者は、同じ薬が予定されていたとしても記録上は対照群に入ります。そのため治療群は、最初の数か月を生き延びた人だけの集まりになります。
素朴な分け方が使っている治療変数を式で書くと、次のようになります。
$$A^{*} = I\left(T_R \le \min(T, C)\right)$$
$T_R$は治療を始めた時刻、$T$は死亡時刻、$C$は打ち切り時刻で、$I(\cdot)$は括弧の中が成り立つ時に1を返す関数である。この変数は、追跡が終わるまでの情報を使って追跡開始時点の群を決めています。未来を見て過去の群を割り当てているわけで、ランダム化試験では起こりえない操作です。
治療開始までの期間は、定義上イベントが起きない期間になります。この期間を治療群の追跡時間に算入すると治療群の分母だけが水増しされ、ハザードは小さく見積もられます。Suissaはこれを不死時間バイアスと呼び、呼吸器疾患や心血管疾患の薬剤疫学研究で同じ誤りが繰り返されてきたと報告しています。今回の架空データでは、治療群に入った1138人が、この構造だけでハザード比を0.42まで押し下げていました。
実際のレセプトデータでは、この誤りは意図せず混入します。抽出条件を「観察期間中に該当薬剤の処方コードが1件以上ある患者」と書いた時点で、既に未来の情報を使った群分けになっているからです。抽出の条件文が、そのまま解析上の仮定になっている状態と言えます。
標的とする試験を先に書き出す
ランダム化試験では、適格性の判定、群の割り付け、追跡の開始が同じ日に起こります。この3つが1点で重なることが、time zeroという言葉の中身です。そして観察データの解析が壊れる時、多くの場合はこの3つが別々の日に分かれています。
HernánとRobinsは、データベースを開く前に、実際には行えない仮想的なランダム化試験の計画書を書くことを勧めています。適格基準、比べる治療戦略、割り付けの方法、追跡の開始と終了、アウトカム、知りたい因果対比、そして解析計画。この7項目を先に書き、それから観察データで各項目をどう再現するかを決めます。
手順としては遠回りに見えますが、time zeroの置き場所が自動的に決まるという利点があります。仮想的な試験で「登録した日に薬を始める」と書いたのであれば、観察データでも薬を始めた日を追跡開始日にするほかありません。設計の段階で決まっていれば、解析の段階で迷う余地がなくなります。私が受託した解析でも、time zeroの定義が文章で残っていない研究は、後から複数の解釈が可能になり、査読で必ず問われる箇所になりました。
7項目のうち、観察データで再現しにくいのは割り付けの方法です。ランダム化に相当するものは観察データに存在しないため、測定された背景で条件付けたうえで割り付けが無作為だったとみなす、という仮定を置くことになります。この仮定が置けるかどうかは、time zeroの揃え方とは別の問題として書き分けます。
開始前の期間を非治療として扱う
起点を直す前に、まず追跡時間の割り当てだけを直してみます。同じ患者を、治療開始前は非治療、開始後は治療として2行に分ける持ち方です。
$$\lambda(t) = \lambda_0(t)\exp\{\beta A(t) + \gamma Z\}, \quad A(t) = I(t \ge T_R)$$
$A(t)$は時刻$t$で治療中なら1、まだ始めていなければ0を返す変数、$Z$は重症度、$\lambda_0(t)$は基準ハザードである。同じ患者が、開始前は非治療の行として、開始後は治療の行として、それぞれの時点のリスク集合に入ります。
この持ち方に直すと、ハザード比は0.42から0.79（95%信頼区間0.72〜0.88）へ動きました。真の値0.70に近づいたものの、区間はわずかに0.70の上側にあります。1つのデータセットで起きたことなので、推定方法の偏りとしては読めません。ここで読み取るべきは0.42と0.79の差の大きさであって、小数第2位の位置ではありません。
開始した日を追跡の起点に置き直す
時変共変量は追跡時間の割り当てを直しますが、比べる集団の作り方までは直しません。重症な患者ほど早く薬を始めるという構造が残っていれば、治療の効果と開始理由の効果は絡んだままです。
Rayが提唱した新規使用者デザインは、薬を始めた日をtime zeroに置き、その時点でまだ始めていない患者の中から背景の近い人を選んで対応づけます。ランダム化試験で全員が同じ日に治療を始めるのと同じ状態を、観察データの中に作ることになります。過去に既に使っていた患者を除く点も重要で、使い続けている人だけを見ると、副作用で早期に中止した患者が最初から抜けた集団を見ることになります。
今回のデータでは606組が対応づけられ、ハザード比は0.83（0.71〜0.96）でした。時変共変量の解析と結果は近く、両方とも0.42とは大きく異なります。傾向スコアマッチングで扱った対応づけの手順自体は同じですが、対応づける時点を治療開始日に揃えるところが新規使用者デザインの中心です。
逐次的に擬似試験を作る
time zeroを1つの時点に固定すると、その前後に薬を始めた患者が解析から消えます。データベース研究では開始時期が数年にわたって散らばるため、この損失は小さくありません。
Danaeiらがスタチンの解析で使った方法は、複数の時点それぞれで擬似的な試験を作り、あとで統合するというものです。今回は3、6、9、12、18、24か月の6時点で、その時期に開始した患者と、その時点でまだ開始していない患者を重症度で対応づけました。
ペア数は3か月時点の170組から24か月時点の48組まで減ります。遅い時点になるほど、まだ薬を始めていない候補者が残っていないからです。統合する時は擬似試験を層として層別Coxに入れ、同じ患者が複数の擬似試験に登場しうるため、分散はクラスタに頑健な推定量で求めます。この点を忘れると信頼区間が狭くなりすぎます。
擬似試験を層に入れるか、1つのデータにまとめて時点を共変量に入れるかは、基準ハザードを共通と考えるかどうかで決まります。開始時期によって背景の患者層が変わっているのであれば、層に入れる方が仮定は弱くなります。
なお、ランドマーク時点で群を決めるランドマーク解析と発想は似ていますが、扱う対象が異なります。ランドマーク解析は治療反応のような追跡中の状態で群を決め、こちらは治療開始そのものを起点に置きます。
Rで確認する
survivalパッケージだけで確認できます。まず素朴な分け方から始めます。
library(survival)
# patients: id, severity, tstart（治療開始時刻、開始しなければNA）, tm, st
# 追跡中に一度でも処方があれば治療群にする持ち方
patients$ever <- as.numeric(!is.na(patients$tstart) & patients$tstart <= patients$tm)
coxph(Surv(tm, st) ~ ever + severity, data = patients)#> coef exp(coef) se(coef) z p
#> ever -0.878 0.416 0.0489 -17.96 <2e-16
#> 治療群 1138人、真の効果は0.70次に、治療開始時刻で行を分割します。tmergeを使うと開始前後の2行に分かれます。
tv <- tmerge(patients, patients, id = id, death = event(tm, st))
tv <- tmerge(tv, patients, id = id, treated = tdc(tstart))
coxph(Surv(tstart, tstop, death) ~ treated + severity, data = tv)#> coef exp(coef) se(coef) z p
#> treated -0.232 0.793 0.0509 -4.56 5.1e-06
#> 6138行（5000人が治療開始で分割されたもの）読み方で迷いやすいのは行数です。6138行は6138人ではなく、5000人のうち1138人が2行に分かれた結果です。人数が増えたわけではないため、標準誤差もこの分割で小さくなるわけではありません。tdcは指定した時刻以降を1にする関数で、開始時刻がNAの患者は全期間0のままになります。
対応づけたデータでは、擬似試験を層に入れて解析します。
coxph(Surv(tm, st) ~ treated + severity + strata(trial), data = matched)
#> treated exp(coef) 0.825 95%CI 0.708-0.960 606組論文にどう書くか
標的試験の7項目を表にして、観察データでの再現方法と並べて載せる書き方が定着しつつあります。とくにtime zeroの定義と、開始日の前後にどれだけの猶予期間を認めたかは、数字が動く箇所なので本文に書いておきます。猶予期間を30日にするか90日にするかで結果がどう変わるかは、感度分析として並べておくと査読での往復が減ります。
除外した患者数も、時点ごとに書きます。逐次的に擬似試験を作った場合は、各時点で何人が適格で何組が対応づいたかを表にすると、図3のような減り方が読者に伝わります。
ここまで読んで、自分のデータベース研究ではどの時点を起点にすべきか判断に迷う方もいるかもしれません。実際には、薬の使われ方、記録の粒度、猶予期間の置き方まで含めて決めることになります。Dr.データサイエンスでは、こうした設計段階のご相談を承っています。
次に確かめること
治療が1回きりの開始ではなく、途中で中止と再開を繰り返す場合は、対応づけだけでは足りなくなります。その時は周辺構造モデルで扱った重み付けと組み合わせることになります。また、time zeroを揃えても測定されていない交絡が残る点は変わらないため、因果推論の前提そのものは別に確かめる必要があります。そして、開始する時期そのものを比べたい場合は、動的治療戦略で閾値ごとの比べ方を扱っています。
参考文献
1. Hernán MA, Robins JM. Using big data to emulate a target trial when a randomized trial is not available. American Journal of Epidemiology. 2016;183(8):758-764.
2. Hernán MA, Sauer BC, Hernández-Díaz S, Platt R, Shrier I. Specifying a target trial prevents immortal time bias and other self-inflicted injuries in observational analyses. Journal of Clinical Epidemiology. 2016;79:70-75.
3. Suissa S. Immortal time bias in pharmacoepidemiology. American Journal of Epidemiology. 2008;167(4):492-499.
4. Ray WA. Evaluating medication effects outside of clinical trials: new-user designs. American Journal of Epidemiology. 2003;158(9):915-920.
5. Danaei G, Rodríguez LA, Cantero OF, Logan R, Hernán MA. Observational data for comparative effectiveness research: an emulation of randomised trials of statins and primary prevention of coronary heart disease. Statistical Methods in Medical Research. 2013;22(1):70-96.
