抗凝固薬をいつ始めるか、透析をどの時点で導入するか、免疫抑制薬をどこまで待つか。臨床で実際に議論されるのは、使うか使わないかではなく、どの状態になったら始めるかです。ところが観察データの解析では、治療あり群と治療なし群の二択に落とし込んでしまうことが多くあります。
周辺構造モデルでは治療を続けるか続けないかという静的な方針を比べました。ここで扱うのは、患者の状態に応じて開始時期が変わる方針です。
比べる対象を規則そのものにする
動的治療戦略とは、その時点までに観測された情報から次の行動を決める規則のことです。「検査値$L$が初めて$\theta$以上になった受診で治療を開始し、以後は継続する」という規則を$d_\theta$と書くと、次のようになります。
$$d_\theta:\quad A_k = I\left(\max_{j \le k} L_j \ge \theta\right)$$
$A_k$は$k$回目の受診時点で治療中なら1を取る変数、$L_j$は$j$回目の検査値、$I(\cdot)$は括弧の中が成り立つ時に1を返す関数です。$\max$を使っているのは、一度でも閾値を超えたら以後は治療を続けるという意味で、規則が過去の情報だけで決まっている点が重要です。
山ができる理由は、治療の効き方が状態によって変わるところにあります。検査値が高い患者では治療が死亡リスクを下げますが、まだ低い患者では、副作用や合併症の分だけ差し引きで不利になります。同時に、治療は検査値そのものを下げるので、早く始めるほど以後の悪化を抑えられます。開始が早すぎれば不利益を先取りし、遅すぎれば悪化を止め損ねるという構造で、その釣り合いが取れる位置に最良の閾値が現れます。
$\theta$を動かせば規則の族ができます。極端な例では、$\theta$を十分小さくすれば最初から全員治療、十分大きくすれば一度も治療しない方針になり、二択はこの族の両端にあたり、臨床的に知りたいのは、どの$\theta$で生存が最もよくなるかどうかという点です。
1人の患者を複数の規則に振り分ける
ここで困るのは、観察データの中に「規則$d_\theta$に従った群」が最初から存在しない点です。患者は医師の判断で治療を始めており、その判断が$d_{0.4}$に一致するか$d_{1.2}$に一致するかは、後から追いかけて初めて分かります。
そこで、追跡開始の時点では全員が全ての規則と矛盾していないと考え、患者を規則の数だけ複製します。そして複製ごとに追跡し、その規則が求める行動と実際の行動が食い違った受診で、その複製を打ち切ります。この手順はターゲット試験エミュレーションで扱った考え方の延長にあります。
同じ1人が、ある規則では最後まで残り、別の規則では途中で消えます。複製同士は独立ではないので、信頼区間は患者を単位とするブートストラップで求めます。
打ち切りが予後と結びついている
この打ち切りは、脱落や転院のような偶発的なものとは性質が異なります。検査値が悪化した患者ほど医師は治療を早く始めるため、閾値の高い規則からは、状態の悪い患者ばかりが先に打ち切られていきます。残った複製は予後のよい患者に偏ります。
そのため、打ち切りの確率で重みを付けます。各受診で規則どおりの行動が取られる確率をモデル化し、その逆数を掛けていきます。
$$W_{ik} = \prod_{j=1}^{\min(k,\,s_i)} \frac{1}{\Pr\left(A_{ij} = d_\theta(\bar{L}_{ij}) \mid \bar{L}_{ij}, Z_i\right)}$$
$\bar{L}_{ij}$は$j$回目までの検査値の履歴、$Z_i$は登録時の重症度、$s_i$は治療を開始した受診の番号です。積を$s_i$で止めているのは、開始したあとは判断の余地がなく、打ち切られる可能性もなくなるからです。重みは各時点までの累積値を使い、最終的な重みを追跡期間全体に当てはめてはいけません。
実務では、閾値を超えた当日に開始する患者はまれです。次回受診までに始めれば規則に従ったとみなす、といった猶予期間を規則の定義に含めておくと、打ち切られる複製が減り、重みも落ち着きます。猶予期間を何か月に置いたかは推定値を動かすので、規則の定義そのものとして書き残します。
重みが大きくなりやすい時は、分子に治療開始の周辺確率を置いた安定化重みを使います。ただし安定化で小さくなるのは重みの散らばりであって、後で触れる正値性の問題そのものは残ります。
重みを外すと順序が逆になる
ランダムで生成した6000例のサンプルデータで検証したところ、重みを付けた推定値は、閾値0.4で0.644、0.8で0.666、1.2で0.662となり、真の値の0.645、0.654、0.644をよく再現していました。最良の閾値も0.8のまま変わりません。
一方、複製と打ち切りだけで重みを付けないと、推定値は0.66から0.69の狭い範囲に固まり、山が消えます。とくに一度も治療しない方針が0.684となり、8つの方針の中でほぼ最良に見えます。真の値は0.483で最下位ですから、順序が逆転しています。予後のよい患者ほど治療せずに済んでおり、その偏りがそのまま方針の成績として読まれた結果です。
二択との差も見ておきます。最初から全員に治療する方針の真の値は0.567、推定値は0.570でした。閾値0.8の規則は真の値で0.654ですから、静的な方針の中で良い方と比べても0.087の開きがあります。100人あたり9人近い差が、開始時期を状態に合わせるかどうかだけで生じている計算になります。二択で比べて0.567と0.483のどちらが良いかを論じても、この差には手が届きません。
極端な規則では推定が壊れる
重みを付けても、一度も治療しない方針の推定値は0.632で、真の値0.483から大きく外れました。原因は正値性です。検査値が高くなった患者が8回の受診を通じて一度も治療されない確率はほとんどゼロで、そのわずかな患者に3171という重みが乗ります。有効標本サイズは45まで落ちており、6000例のデータを使いながら実質45例で推定していることになります。
私なら、有効標本サイズが元の1割を切った規則は結果として報告しないことを検討します。信頼区間を付ければ誠実に見えますが、区間そのものが少数例の偶然に左右されており、幅を広げても正値性の破れは埋め合わせられないからです。報告するのは、データが実際に支えている閾値の範囲だけに限ります。今回であれば0.0から1.2までです。
重みを99パーセント点で打ち切る対処もよく使われます。分散は下がりますが、代わりに偏りが入るので、打ち切りは推定の不安定さを隠す道具にはなりません。有効標本サイズが45まで落ちている状況では、打ち切っても値が信用できるようにはならず、その規則を比較対象から外す方が適切と判断されます。
Rで確認する
まず、各受診で治療を開始する確率のモデルを当てはめます。未治療のあいだの受診だけを行にします。
dd <- do.call(rbind, lapply(seq_along(rows), function(i) {
r <- rows[[i]]; s <- which(r$A == 1)[1]
k <- if (is.na(s)) seq_len(r$last) else seq_len(s)
data.frame(L = r$L[k], sev = r$sev, A = r$A[k]) }))
start_model <- glm(A ~ L + sev, family = binomial, data = dd)
coef(start_model)#> (Intercept) L sev
#> -1.53 1.64 0.40次に、規則が求める行動を作り、逸脱した受診で打ち切ったうえで、各時点までの累積重みを持たせます。
prescribe <- function(L, theta) as.numeric(cummax(L) >= theta)
expand <- function(r, theta) {
pres <- prescribe(r$L, theta)
dev <- which(r$A != pres)[1]
last <- if (is.na(dev)) r$last else dev - 1
if (last < 1) return(NULL)
s <- which(r$A == 1)[1]; if (is.na(s)) s <- r$last
w <- cumprod(1 / ifelse(pres == 1, r$p, 1 - r$p))[pmin(seq_len(last), s)]
data.frame(k = seq_len(last), w = w,
d = c(rep(0, last - 1), if (is.na(dev)) r$dead else 0))
}あとは受診ごとに加重した死亡割合を求め、掛け合わせます。
pp <- do.call(rbind, lapply(rows, expand, theta = 0.8))
haz <- tapply(pp$w * pp$d, pp$k, sum) / tapply(pp$w, pp$k, sum)
prod(1 - haz)
#> 0.666 （閾値2.0では0.581、一度も治療しない方針では0.632）読み違えやすいのは$\theta$を大きくした時の挙動です。閾値を上げるほど治療開始が遅くなるので推定値は単調に下がると思いがちですが、実際には山ができます。開始が早すぎる方針は、まだ検査値の低い患者にも治療の不利益を与えているためです。cummaxを使っている点も注意が必要で、これを外すと検査値が一時的に下がった受診で治療を中止する規則になり、比べている対象そのものが変わります。
論文にどう書くか
比べた規則の一覧を、閾値の値とともに表にします。規則ごとに、複製の数、打ち切られた数、有効標本サイズ、重みの最大値と分布の要約を並べると、どの規則がデータに支えられているかが読者に伝わります。そのため、重みを打ち切った場合は、その閾値と、打ち切りの前後で推定値がどう動いたかも書きます。
また、図2のように山の頂上付近は平らで、0.8と1.2の差は真の値でも0.010しかありません。選んだ1点ではなく曲線の形を見せる方が、読者の判断材料になるため、最良の閾値を1つ選んで報告する書き方には注意が必要です。
ここまで読んで、自分の研究でどこまで細かい閾値を比べられるか判断に迷う方もいるかもしれません。実際には、受診間隔、検査の欠測、開始の猶予期間の置き方まで含めて決めることになります。Dr.データサイエンスでは、こうした解析設計のご相談を承っています。
次に確かめること
今回は測定された重症度と検査値で交絡が全て説明できると仮定しました。この仮定そのものは因果推論で扱ったとおり、データからは確かめられません。また、検査値の履歴をどう持つかで規則の定義そのものが変わるため、時変共変量のデータの持ち方も併せて確認しておくと迷いが減ります。
参考文献
1. Robins JM. A new approach to causal inference in mortality studies with a sustained exposure period. Mathematical Modelling. 1986;7(9-12):1393-1512.
2. Hernán MA, Lanoy E, Costagliola D, Robins JM. Comparison of dynamic treatment regimes via inverse probability weighting. Basic & Clinical Pharmacology & Toxicology. 2006;98(3):237-242.
3. Cain LE, Robins JM, Lanoy E, Logan R, Costagliola D, Hernán MA. When to start treatment? A systematic approach to the comparison of dynamic regimes using observational data. International Journal of Biostatistics. 2010;6(2):Article 18.
4. Murphy SA. Optimal dynamic treatment regimes. Journal of the Royal Statistical Society Series B. 2003;65(2):331-355.
5. Young JG, Cain LE, Robins JM, O’Reilly EJ, Hernán MA. Comparative effectiveness of dynamic treatment regimes: an application of the parametric g-formula. Statistics in Biosciences. 2011;3(1):119-143.
